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Continuing with the example of a transversal filter, suppose we are given a set 
of observations represented by the tap-input vector u(i), where the time index i = 
1 , 2, . . . , n. We assume that the tap-weight vector w of the transversal filter is held 
constant for the entire observation interval. At time f, we define the estimation error 
as 

e(i) =d(i) - w*u(f) (1.20) 

where d(i) is the corresponding value of the desired response. In this case, we may 
define the cost function as 

^(w,«) = i| e (/)p 

T (1.21) 

= 2k(/)-wu(oi 2 

Note that whereas the ensemble-averaged cost function J(w) in Eq. (1.19) depends 
only on the tap-weight vector w, the time-averaged cost fiinction ^(w, n) depends 
on w and the observation interval n. Accordingly, the minimization of ^(w, n) with 
respect to w yields a solution for the tap-weight vector that varies with the observa- 
tion interval. This second approach forms the basis of the method of least squares. 

The cost functions 7(w) and ^(w, n) are both convex with a unique minimum 
point. Accordingly, their use yields a unique solution for the tap-weight vector of the 
transversal filter. A qualification in the context of the method of least squares, how- 
ever, is in order. In particular, in Eq. (1.21) it is assumed that the number of obser- 
vations n is larger than the number of tap weights constituting the vector w; that is, 
we have an overdetermined system with more equations than unknowns. 

A limitation of second-order statistics (e.g., the mean-square-error criterion) 
is that they are phase blind. We may overcome this limitation by the use of a non- 
linear cost function. By so doing, the filter is enabled to extract information 
(particularly phase information) from the input signal in a more efficient manner. 
For this to be possible, however, the input signal must have non-Gaussian statistics. 
The use of such an approach provides the basis for an important class of nonlinear 
adaptive filtering algorithms that can perform blind deconvolution, blind in the sense 
that the algorithms do not require a desired response. The issue of blind deconvolu- 
tion (and its use in blind equalization) are covered in Chapter 20. 



TL7 APPLICATIONS 

The desirable features of an adaptive filter, namely, the ability to operate satisfacto- 
rily in an unknown environment and also track time variations of input statistics, 
make the adaptive filter a powerful device for signal-processing and control applica- 
tions. 

Indeed, adaptive filtering has been successfully applied in such diverse fields as 
communications, radar, sonar, seismology, and biomedical engineering. Although 
these applications are indeed quite different in nature, nevertheless, they have one 
basic common feature: An input vector and a desired response are used to compute 
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an estimation error, which is in turn used to control the values of a set of adjustable 
filter coefficients. The adjustable coefficients may take the form of tap weights, 
reflection coefficients, or rotation parameters, depending on the filter structure em- 
ployed. However, the essential difference between the various applications of adap- 
tive filtering arises in the manner in which the desired response is extracted. In this 
context, we may distinguish four basic classes of adaptive filtering applications, as 
depicted in Fig. 1.7. For convenience of presentation, the following notations are 
used in this figure: 

u = input applied to the adaptive filter 
y = output of the adaptive filter 
d — desired response 
e — d — y = estimation error. 

The functions of the four basic classes of adaptive filtering applications depicted 
herein are as follows: 

L Identification [Fig. 1.7(a)], The notion of a mathematical model is fundamental 
to sciences and engineering. In the class of applications dealing with identi- 
fication, an adaptive filter is used to provide a linear model that represents the 
best fit (in some sense) to an unknown plant. The plant and the adaptive filter 
are driven by the same input. The plant output supplies the desired response for 
the adaptive filter. If the plant is dynamic in nature, the model will be time 
varying. 

II. Inverse modeling [Fig. 1.7(b)]. In this second class of applications, the function 
of the adaptive filter is to provide an inverse model that represents the best fit 
(in some sense) to an unknown noisy plant. Ideally, the inverse model has a 
transfer function equal to the reciprocal (inverse) of the plant's transfer func- 
tion. A delayed version of the plant (system) input constitutes the desired 
response for the adaptive filter. In some applications, the plant input is used 
without delay as the desired response. 

ID. Prediction [Fig. 1.7(c)]. Here the function of the adaptive filter is to provide the 
best prediction (in some sense) of the present value of a random signal. The 
present value of the signal thus serves the purpose of a desired response for the 
adaptive filter. Past values of the signal supply the input applied to the adaptive 
filter. Depending on the application of interest, the adaptive filter output or the 
estimation (prediction) error may serve as the system output. In the first case, 
the system operates as a predictor; in the latter case, it operates as a prediction- 
error filter. 

IV. Interference cancelling [Fig. 1.7(d)]. In this final class of applications, the 
adaptive filter is used to cancel unknown interference contained (alongside an 
information-bearing signal component) in a primary signal, with the cancella- 
tion being optimized in some sense. The primary signal serves as the desired 
response for the adaptive filter. A reference (auxiliary) signal is employed as the 
input to the adaptive filter. The reference signal is derived from a sensor or set 
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Figure 1.7 Four basic classes of adaptive filtering applications: (a) class I: 
identification; (b) class H; inverse modeling; (c) class ffl: prediction; (d) class IV: 
interference cancelling. 
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of sensors located in relation to the sensor(s) supplying the primary signal in 
such a way that the information-bearing signal component is weak or essentially 
undetectable. 



In Table 1 . 1 we have listed some applications that are illustrative of the four 
basic classes of adaptive filtering applications. These applications, totaling twelve, 
are drawn from the fields of control systems, seismology, electrocardiography, com- 
munications, and radar. 6 They are described individually in the remainder of this 
section. 



TABLE t.1 APPLICATIONS OF ADAPTIVE FILTERING 


Class of adaptive filtering 


Application 


I. Identification 


System identification 




Layered earth modeling 


H. Inverse modeling 


Predictive deconvolution 




Adaptive equalization 


m. Prediction 


Linear predictive coding 




Adaptive differential pulse-code modulation 




Autoregressive spectrum analysis 




Signal detection 


IV. Interference cancelling 


Adaptive noise cancelling 




Echo cancellation 




Radar polarimetry 




Adaptive beamforming 



System Identification 

Systemidentification is the experimental approach to the modeling of a process or a 
plant (Astrom and Wittenmark, 1990; Soderstrom and Stoica, 1988; Ljung, 1987; 
Ljung and Soderstrom, 1983; Goodwin and Payne, (1977). It involves the following 
steps: experimental planning, the selection of a model structure, parameter estima- 
tion, and model validation. The procedure of system identification, as pursued in 
practice, is iterative in nature in that we may have to go back and forth between 
these steps until a satisfactory model is built. Here we discuss briefly the idea of 
adaptive filtering algorithms for estimating the parameters of an unknown plant 
modeled as a transversal filter. 

Suppose we have an unknown dynamic plant that is linear and time varying. 
The plant is characterized by a real-valued set of discrete-time measurements that 
describe the variation of the plant output in response to a known stationary input. 
The requirement is to develop an on-line transversal filter model for this plant, as il- 
lustrated in Fig. 1.8. Hie model consists of a finite number of unit-delay elements 
and a corresponding set of adjustable parameters (tap weights). 

Let the available input signal at time n be denoted by the set of samples: w(n), 
u(n — 1), . . . , u(n — M + 1), where M is the number of parameters in the 

6 For additional applications of adaptive filtering, see Widrow and Stearns (1985), Cowan and 
Grant (1985), and the special issues on adaptive filters listed in the References and Bibliography at the 
end of the book. 
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Comparison of RLS, LMS, and Sign Algorithms 
Tracking Randomly Time-Varying Channels 
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&£strach ~ The performance of adaptive FIR filters governed 
by the recursive least-squares (RLS) algorithm, the least mean 
square (LMS) algorithm, and the sign algorithm (SA), are com- 
pared when the optimal filtering vector is randomly time-varying . 
The comparison is done in terms of the steady-state ^excess 
mean-square estimation error £ and the steady-state mean-square 
weight deviation, i;. It is sfioTvn that £ does not depend on the 
spread of eigenvalues of the input covariance matrix* in the 
cases of the LMS algorithm and the SA, while it does in the 
case of the RLS algorithm. In the three algorithms, q is found 
to be increasing with the eigenvalue spread. The value of the 
adaptation parameter that minimizes £ is different from the one 
that minimizes //. It is shown that the minimum values of £ and rj 
attained by the RLS algorithm are equal to the ones attained by 
the LMS algorithm in any one of the three following cases; (1) if 
f~£? has equal eigenvalues, (2) if the fluctuations of the individual 
-elements of the optimal vector are mutually uncorrelated and 
"Aave the same mean-square value, or (3) if E is diagonal and the 
- fluctuations of the individual elements of the optimal vector have 
Mhe same mean-square value. Conditions that make the values of 
and jj of the LMS algorithm smaller (or greater) than the ones 
%jof the RLS algorithm are derived. For Gaussian input data, the 
^minimum values of £ and ?/ attained by the SA are found to exceed 
-^thc ones attained by the LMS algorithm by 1 dB independently 
~r=of /? and the mutual correlation between the elements of the 
s optimal vector. 

%1 I. Introduction 

1HE performance of an adaptive filter depends mainly 
on the algorithm used for updating the filter weights. 



time-varying channel. While a reasonable progress fuis 
achieved in the tracking analysis of the LMS algorithm and the 
SA (e.g., [6], [7], [9], [12]), the available tracking analysis of 
the RLS algorithm [l]-[3], [14] is far from being complete. In 
[1], the RLS algorithm is studied when the signal characteristic 
has a step change, and it is shown that the algorithm converges 
exponentially with a time constant equal to one half of 
the reciprocal of the adaptation parameter. In [2], the RLS 
algorithm is studied in the case when the fluctuations of 
the individual channel parameters are Gaussian, mutually 
uncorrelated, and of the same mean-square value. Under these 
assumptions, it is shown [2, p. 1 102] that the RLS algorithm 
has the same tracking capability as the LMS algorithm when 
the input covariance matrix has equal eigenvalues. In [3], 
it is shown that the performance of the RLS algorithm is 
worse than that of the LMS algorithm in. the case of tracking 
a deterministic chirped sinusoid buried in additive white 
noise. 

f The j >rescni paper provides a tracking analysis of t he RLS 
alg orithm that do es not include assumptions on the typ e 
of distribution and mutual correlation amon^j he chann el 
parameter s. The paper compares the RLS, LMS, and SA 
algorithms from the points of view of the steady-state exce ss 
mean-square estimation error £ a nd the steady- state mea n- 
square weight deviation, rj r As will be shown in the paper, 
the value of the adaptation parameter /x, that minimizes £, 



OThree well-known algorithms are: the recursive least-squares f is generally d^rent from the one that minimizes tj. Thus, 
H(RLS) algorithm [l]-[5], the least mean-square (LMS) algo-J lhe optimum value of ji depends on whether £ or r, is 
rithm (6], [I0H12], and the sign algorithm (SA) [7H°]> 35 ^ nm ^ Pe rformance index * P^ 1 * denves 
"Convergence analyses of the three algorithms in the case of expressions of i\ for the three algorithms and generalizes 
stationary adaptive filtering lead to the following conclusions. 



With appropriate choice of the adaptation parameter, the RLS 
and LMS algorithms are exponentially convergent [5], [10] 
whereas lhe SA is linearly convergent [8]. This makes the 
convergence of the SA the slowest when the initial weight 
setting is far from the optimum. The convergence rate of 
t he LMS algorithm decreases as the spread of eigenvalue s 
of th e input covariance matrix increases [II] whereas t he 
convergence rate of the RLS algorithm is independent of the 
eigenvalu e sprea d [5]. 

The purpose of this paper is to give a quantitative com- 
parison of the three algorithms when tracking a randomly 

Manuscript received March 19, 1993; revised March 1, 1994. The associate 
editor coordinating the review of this paper and approving it for publication 
was Dr. Fuyun Ling. 

The author is with the Military Technical College, Kobry El-Kobba, Cairo, 
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IEEE Log Number 9404766. 



previously derived expressions of £. For each algorithm, the 
value of /i that minimizes 77 and the one that minimizes £ 
are derived. The paper compares the three algorithms on the 
base of the minimum values of 77 and £ attained by each 
algorithm. Conditions that make the tracking properties of 
the LMS algorithm better (or worse) than those of the RLS 
algorithm are derived. The paper is organized as follows. 
Section II is concerned with the problem formulation and 
the assumptions used throughout the paper. Sections III-V 
are concerned with the derivatio n of £ a nd^ of t he RLS 
algorithm, th e LMS algorithm ;"and the sign algorithm, re- 
spectively. Section VI compares the minimum values of £ 
and t; of the three algorithms and provides a decision tree 
for the choice between the RLS and LMS algorithms. Section 
VII provides computer simulations that support the derived 
theoretical results. The conclusions of the paper are given in 
Section VIII. 
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\IL/ Problem Formulation 

Let kk,£k- an d a k denote the weight vector of the adaptive 
filter at discrete time k, the observation vector, and the desired 
filter output, respectively. The estimation error e k is given by 

e* = a* - fcjjTak (2J) 

where k[ is the transpose of h k . The RLS algorithm is given 
by [IH5] 

=4fc + Affit+iacfc, 0 < < 1 (2.2) 



1 



(2.3) 



In (23), Pj is an arbitrary symmetric positive definite matrix. 
The quantity (1 - p) is usually referred to as the forgetting 
factor of the RLS algorithm. The LMS algorithm is given by 
[6] 



: h k + M£* e fc, 0 < n < iiQ 



(2.4) 



where [mq is a positive number that depends on the statistics 
of x k . Finally, the SA is given by [7]-[9] 

= hk + V£k sgn ( e *) » M > 0. (2.5) 

|r ~* Let the relation between the desired filter output a* and the 
observation vector, be modeled by 

<ik = c^£ fe + ftft * (2.6) 

where c fc and are the model parameter vector and the model 
noise, respectively. Examples of adaptive filtering applications 
that follow the model (2.6) are adaptive system identification 
and adaptive echo canceling [15J. The vector c fe is assumed 
randomly time-varying. Let the increment of Cj. be denoted by 



(2.7) 



The following assumptions are used throughout the paper; 

^ A.1) The sequences {x k },{b k }^ and {d k } are mutually 
independent. 

^ A.2) {dk} is a stationary sequence of independent zero- 
mean vectors, and Q = B(d k d[) is the covariance 
matrix of the increments. 
A.3) The sequence {x k } is stationary, zeromean, and the 
covariance matrix R = E(x k x^[) is positive definite. 
A.4) {b k } is a stationary sequence of independent ze- 



romean random variables, with variance of. 



These assumptions will be used to analyze the tracking of 
the RLS and LMS algorithms. The tracking analysis of the 
SA will use the assumptions (A.1HA.4) along with the 
assumption that the sequences {b k } and {x k } are Gaussian. 
These assumptions are the same as the ones considered in [7], 
[9], and [13]. Another model of variation of c,, that has been 
considered in the context of nonstationary adaptive filtering is 
a Markovian model f61 that can be expressed as 



model is close to the model expressed by the assumption (A.2). 
Indeed, it is easy to show that this model satisfies 

ml^l/E^d,) = (1 - a)/2 » 0 

which means that (J& and are nearly uncorrected. It 
should be noticed that the assumptions (A.1HA.4) do not 
imply that the individual elements of x k have the same mean- 
square value. Therefore, (A.1HA.4) allow the case when the 
individual elements of x k are derived from different signals 
as well as the case when x k is made of samples of the same 
signal. 

The steady-state excess mean-square estimation error £ is 
defined by 

VmE(el)-crl 

k — *oc 

Let the instantaneous weight deviation from the optimal be 
denoted by 




h k -c k . <^ (2.9) 

The steady-state mean-square weight deviation ?/ is defined by 



\im E{\\v k f) 



where j|v fc j| 2 = v^v k is the squared norm. In the following 
sections, we derive expressions of £ and 77 for each algorithm. 

f^m^ Tracking Properties of the RLS Algorithm 
Denote 



Ms 1 . 



(3.1) 



Equations (2.3), (3.1), and the matrix i nversion lemma imply 
(i'Bk+i = (l-/0& + MM-f. J (3-2) 

t _ 

From (2.1), (2.2), (2.6), (2.7), (2.9), and (3.1) one obtains 

Sk+i = + V&iliXkibk - 3£fr£*) " dk- (3.3) 

Premultiplying both sides of (3.3) by R^+i and then using 
(3.2), one obtains 

Due to (3.2), the sequence {R k } depends only on the sequence 
Then (A.l) implies that {R k } is independent of the 
sequences {b k } and {d^}. Due to (3.3), v k is a function of 
{bj,dj,xj);j < k}. Therefore, (A.l), (A.2), and (A.4) imply 
that b k and are independent of v k . Squaring (3,4), using the 
independence of b k and d k on v^x^R^^ and Bx+i* using the 
fact that bk and d^ are zeromean and that B{b\) — a\ and 
E{&&) = Q imply that 

= (3.5) 
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Therefore, for small values of we can use the approximation 
that 7^. = R in the steady state. In such a case, for large values 
of k, (3.5) can be approximated by 

= (1 - tifRE{v^r[)R + + (3.7) 

Premultiplying and postmultiplying the two sides of (3.7) by 
R~~ l yields 

E(v M v^) = (1 - /i) 2 E(a^) + M 2 ^fi" 2 + Q- 0-8) 
Since 0 < /z < 1, then (3.8) implies that 



the value of ji that minimizes £ RLS , and the minimum value 
of f RLS , are, respectively, given by 

■ MS. WO.®} 1 '* 

Nat J 



(3.15) 
(3.16) 



lim ^(jiaZ") = ^i- 



W h IC x +ir l g\. (3.9) 



From (3.11) and (3.15), it is obvious that fif LS is different 
from (i^* In the particular case when R = cJ, /x R LS = m£ ls . 
When Q = cj> R LS = >/ctr(£)/(A^) and - 

ycN/(a% tr(R~ 1 )}. In the latter case, fif LS does not depend 
on the eigenvalue spread of R while /x RLS tends to zero as the 
minimum eigenvalue of R. lends to zero. 



The steady-state mean-square weight deviation is then given 
by 



V™- 5 £ lim EUfaf) 

k — »oo 
1 



rSJ where tr (.) denotes the trace of (.). The first term on the RHS 
of (3.10) is due lo th e plant noise while the second one is 
f7 due to the plant nonstationarity . The first term is increasing 
j ^ in ft whereas the second one is decreasing in jr.O < /* < 1. 
I i Assuming that // <C 1» as is usually the case in applications, 
;2 then the value of /z, denoted by m£ xs , that minimizes tj 1 **^, 
and the minimum value of ?j aLS , denoted by ??^f , are, 
respectively, given by 



tr(Q) 



,1/2 



'/rain 



<T b] /tv(Q)tT(R^y 



r Now, to find the steady-state mean-square estimation error 
^rls^ we $ | ia j| use ^ e f 0 j} ow j n g assumption: ^ 

j^^fr ffi^ ) and :n fe are statistically independent.* v^vi^ 
This assumption has been used by many authors (e.g., [7], [8], 
[13]) and is approximately valid at small values of /i. Due to 
(2.1), (2.6), and (2.9), e k = b h - y£x k . Then (A3) and the 
fact that 6jt is zeromean and independent of z/£ a^. yield 

E(el) = a? + BOdfBfafcsDa). (3.13) 

Equations (2.8), (3.9), and (3.13) imply that the steady-state 
excess mean-square error of the RLS algorithm is given by 



1 



[A^+M-HrfQR)] 



(3.14) 



where A r is the dimension of ^ Hie first term on the RHS 
of (3.14) is due to the plant noise while the second one is 
due to the plant nonstationarity. In the particular case when 
Q — cL, with J being the identity matrix and c being an 
arbitrary positive number, the result (3.14) reduces to the one 
obtained in [2, p. 1 102]. From (3.14), assuming that y, < 1, 



(1v^ Tracking properties of the LMS Algorithm 

Under the assumptions (A.1)-(A.5X and with ti satisfying 
fi tr (R) -C 1, it has been shown in [7, p. 2053] that the steady- 
state excess mean-square error of the LMS algorithm is given 
by 

f LMS = Jfrof tatffl + M- 1 tr(Q)]. (4.1) 

In the particular case when Q — cJ, (4.1) reduces to the result 
obtained in [6]. Due to (4J), the value of \l that minimizes 
£ LMS , and the minimum value of £ LMS , are, respectively, 
given by 

1/2 

(4.2) 



LMS r tr £) 



(3.11) 



(3.12) 



fSi S =^tr(fi)tr(Q). 



(4.3) 



In order that /*£ MS satisfy the condition ^ MS tr (R) < 1, 
one should have 



^tr(£)tr(Q)<<7 6 . 



(4.4) 



This condition is usually called the condition of slow variations 
[12]. 

® In the following, we derive the steady-state mean-square 
weight deviation, t^ us . From (2.1), (2.4), (2.6), (2.7), and 
(2.9) one has 

Uk+i = (Z - f*£k£k)Vh + V>bkX k - (4.5) 

Squaring (4.5) and using the assumptions (A.1MA.5) yield 

■Efefc+iUt+i ) = + i^vlSL + Q (4.6) 

where 

= £[(/ - - 

+ ^(asfJBfefciJ)^). (4-7) 

For small ^, the last term on the right-hand side of (4.7) can 
be neglected with respect to the second term and then (4.6) 
and (4,7) yield 

£(2k+l£*H-l) 
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In the steady state, E(v M ^^) = E(v k y£) and then (4.8) and the corresponding optimum value of \i are* respectively, 
gives given by 



Km [»D + Efa*£)M = wtR + H^Q. (4.9) « tT m tr {Q) 

k~*oo 4 



+ 



|^ 2 ir(7?)tr(Q)-f ^tr 2 (7?)tr 2 (Q) 



1/2 



(5.1) 



It should be noticed thai (4.1) can be obtained by taking the 
trace of both sides of (4.9). Now, since the covariance matrix 

R is symmetric and positive definite, it can be decomposed as ^sa _ ^/tr(Q)/tr (J?). (5.2) 

R = UAU (4.10) Equation j5 2) i s the same as the one derived in [9]. Under 

TT . , - . , the condition (4.4) of slow variations, (5,1) reduces to 
where U_ is an ortbonormal matrix and 

A^diag^,.,^) (4.11) ^ n = ^*ir(R)tr(Q)/2. (5.3) 



It has also been shown in [7, f51 1 that 



where 



with Xw'^n being the eigenvalues of R. Premultiply- 

ing and postmultiplying both sides of (4.9) by if and U, y/2pKp\S^ + M'A" 1 ]/^ ^ ^ + Q'A" 1 (5.4) 
respectively, yields 

where A, Q', and S' are defined by (4.1 1), (4.13), and (4.14), 

AS' + $'A = V>&b A + M Q' (4-12) respectively. Taking the trace of both sides of (5.4) and using 

(4.14), then 

, A T V^M tr (S)/cr 6 = AT M 2 + tr (Q'A -1 ). (5.5) 
$*U T QU (4.13)1 

5' 4 [/ r SLT (4 14)1 Due to (4 * 10) ^ (4 ' 13) ' ^(S'A" 1 ) = tr^fi- 1 ). TTien 

' (5-5) and (4.15) imply that the steady-state mean-square weight 

S - lim BfefcuT). (4.15)j deviation of the SA is given by 

fc— *oo 

Postmultiplying both sides of (4.12) by A -1 and then taking r ? SA = ^Eihkf) 

the trace of both sides yield = J^a&N + ^ tr (QRT 1 )}. (5.6) 

tr(S') = {[tiNtf + fi- 1 tr(g A- 1 )], y (4.16) ■ ^ ^ ya]ue rf ^ ^ mMmizes ^sa and ^ minimum 

Due to (4.10) and (4.13), tr (Q'A" 1 ) = tr (^IT 1 ), and then vaIue of ^ «* respec tively, given b y 

(4.14H4.16) imply that the steady-state mean-square weight SA _ / (np-M/AT j\ 

deviation of the IMS algorithm is given by M * " V tr vi & )/J { ' } 

fc— +00 / . 

2 fc ^ y the eigenvalue spread of R while ;z^ A does. In the particular 

The value of ft that minimizes t? lms and the minimum value case when # = cj\ fif A = wj A - It is worth mentioning that, 

of t/ lms are, respectively, given by due to (4.2), (4.18), (5.2), and (5.7), the optimum vjtfues of 

p, of the LMS algorithm and the SA are related to each other 

1% MS = \/tr(QR- l )/(Naf) (4.18) according to 

t I T~ ,,SA _ „LMS ..SA _ .IMS ( z q\ 

T^=a b y/NtT(QRr 1 ). (4.19) H ~ ° b H ■■ I'v ~ • ^ 

From (4.2) and (4. 1 8), it is seen that $ MS does not depend on » (fo) Comparison of the Tracking 

1 the eigenvalue spread of R while ^ MS does. In the particular J PROPERTIES OF THE ALGORITHMS 

lease when R = cJ,^ MS = f /y AS . )y / Table I summarizes the above derived minimum mean- 

"V*""" ~~ square errors, minimum mean-square weight deviations, and 

V. TRACKING PROPERTIES OF THE SA optimum adaptation parameters of the three algorithms 1 From 

, . . - "is table, one observes that £^ s , /*£ MS , £„t„. and /*f A do not 

The analysts g,ven ,n Sect,on ffl and IV does noMnclude m ^ of eigenvaIues € of & whereas \*£ and 

any assumption on the type of distribution of x h and b k . In rl S ^ „ tK _ u , _ „ ^ f 
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TABLE I 

Minimum Mean-Square Error, Minimum Mean-Square 
Weight Deviation, and Optimum Adaptation Parameters 
of the RLS Algorithm, LMS Algorithm, and SA 





RLS Algorithm 


LMS Algorithm 


Sign Algorithm 










o\ / Ntr(QR) 
o 


<j\ /trCR)tr(Q) 

D - 


<r b Vittr{R)trtQ)/2 


\in 








<r./ trtQHrCR" 1 ) 


<r./ Ntr(QR~* } 
b — 






r tr((JR) ,1/2 






L ^ 1 

0 




/ tr{Q)/tr(R) 


• 


r tr(Q) nl/2 






l ^trtR" 1 ) ' 


J trCQFf 1 }/CN<^) 


/ trCQR -1 )/N 



weight deviation increases as the minimum eigenvalue of R 
decreases. From Table I one has 

& " Ltr 



: and 



RLS 

*min 

''nun 



Cg)*r(g)J 



tr(g)tr(2r ) 



(6.1) 



A r tr(Q£ _1 ) 



1/2 



(6-2) 



From (6.1) and (6.2), ffi = ^ «nd ^ = ^ ^ ^ 

one of the three following cases: (1) if R = cl, with c being 
an arbitrary positive number, (2) if Q = of, or (3) if J? is 
diagonal and the diagonal elements of Q are equal. In these 
three cases, the LMS algorithm is more recommended than the 
RLS algorithm due to the technical simplicity of the former. 

From (6.1), f*j* n s may be greater or smaller than £^ s , 
depending on whether Nii(QR) is greater or smaller than 
tr{Q)tr(ii), respectively. Similarly, from (6.2), r£j£ may 
be greater or smaller than ?7^ s , depending on whether 
tr(Q)tr(^ _1 ) is greater or smaller than N%i(QR~ l ), 
respectively. In the following, we give a geometrical meaning 
of these conditions. Let A 1: A2, * - * , Aiv; Ai > A 2 > ■ • > 
A a- > 0 be the eigenvalues of the matrix R, let u x , t^, ■ * * , u N 
be the associated eigenvectors, respectively, and let 

7i = ^((&) 2 ), * = L2 5 -",A r (6.3) 

denote the mean square of the projection of along the 
direction of u-. In the appendix, we show that 

p NtT{QR)-ti{Q)tr{R) 

5 N-l N 

= EE (^-^)(7i-7i) (6-4) 

I Nlr(QR- l )-tr{Q)tr(ir l ) 

i X—l N 

= EE (^r 1 -A- 1 )(7 i -7i). (6.5) 



Since A a - > \ 5 for i < j, then due to (6.1), (6.2), (6.4), 
and (6.5), a sufficient condition of and ^p n s to be 



greater (resp. smaller) than £j^ s and 77^,, respectively, is 
that 7^ > 7j (resp. 7$ < 7^) for i < j. Thus, a sufficient 
condition for the LMS algorithm tracking to be better (resp. 
worse) than the RLS algorithm tracking is that the mean-square 
projection of the optimal vector increment along the direction 
of an eigenvector of J? be increasing (resp. decreasing) with 
the associated eigenvalue. The result matches the well-known 
fact [11] that the adaptation speed of the LMS algorithm along 
the direction of an eigenvector of R is increasing with the 
associated eigenvalue, while the adaptation speed of the RLS 
algorithm is uniform with respect to the direction [1], [2], [5]. 

The above conditions can be used to choose between the 
RLS and LMS algorithms when the fluctuation covariance 
matrix # is known. Now, we derive a condition that can 
be used when Q is unknown. Since tr(QR) — E(d[Rd k ) 
and tr(Q) = £(&), then A min tr(Q) < tr (QR) < 
A m ax tr (Q) where A m i n and A max , respectively, denote the 
minimum and maximum eigenvalues ofR Then, (6.1) implies 
that 



V^Utl< S < VNX max /tr(R). 



(6.6) 



Similarly, (6.2) implies that 



/ ^RLS i 

VWrCfl- 1 )/^ < ^ < V^maxtr^- 1 )/^. (6.7) 



The lower (resp. upper) bounds in (6.6) and (6.7) are at- 
tained when the parameter fluctuation vector d: k has the same 
direction as the eigenvector associated with the minimum 
(resp. maximum) eigenvalue. Equations (6.6) and (6.7) can 
be used to choose between the RLS and the LMS algorithms 
when the matrix Q is unknown as follows. When f is the 
primary performance index, it is reasonable to choose the 
LMS algorithm when the geometric mean of the lower and 
upper bounds in (6.6) is greater than or equal to 1, i.e., to 
use the LMS algorithm if N^X m i n X mBX > tr(fl) and to 
use the RLS algorithm otherwise. Similarly, when r] is the 
primary performance index, it is recommended to use the LMS 
algorithm if tr {RT 1 ) VA m i n A max > N and to use the RLS 
algorithm otherwise. The above arguments about the choice 
between the RLS and LMS algorithms when Q is known and 
when Q is unknown are summarized in the decision tree given 
in Fig. 1. 

Finally, from Table I, for Gaussian bk and and for all 
values of R and Q, one has 



d-SA ycLMS _ „SA /„LMS _ /TTo 



(6.8) 



Thus, and i/^a exceed ^lms and by 
about 1 dB. Then, for Gaussian 6* and the choice of SA 
versus LMS algorithm depends on how a 1-dB performance 
difference improvement compares to the difference in the 
implementation complexity of the two algorithms. 



VII. Simulation Results 

Hie simulations are done for a case of adaptive identification 
of a randomly rime-varyine plant described bv (2.6) where Ck 
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H - cl 7 or Q - el? 



Is R diagonal? and are 
the~di agonal elements 
of Q equal? 



use LMS 



Mean 
square 



Wiat is the primary 
performance index? 



^NtrtOjp-trtQKrlR) l "° 



r 



RLS 



a 



use use 
LMS RLS 



Mean 
square 
weignt 
deviation 



"M^minNza* * *«■«» 7 



use use 
LMS RLS 



use 
LMS 



Rg. 1. Choice between the RLS and the LMS algorithms, 

is the plant-parameter vector, bk is the plant noise, and Xj. is 
the observation vector given by 



(7,1) 



with Zk being the plant input Hie sequences {bk} and 
are mutually independent, zero mean, and Gaussian. Hie 
sequence {bk} is white. The sequence {x k } is a stationary 
autoregressive sequence given by 



Xk = /ten + \/l-fPvk> 0 < /3 < 1 



(7.2) 



where is a zero mean unity variance white Gaussian 
sequence. The parameter 0 controls the degree of correlation 
of the sequence {xfc}; the greater 0 is, the stronger the 
correlation. TTie used model of variation of the plant-parameter 
vector Cj. is 

Cfr+i = &k + ik> dk = fai.jfe, a>2jh • • ■ , wat,*) 7 * (73) 

where {w 2 ,Ar}.* -** and {wjvjt} are zero mean white 

Gaussian sequences with the same variance a\. The sequences 
{toi t fc} 5 {«/ 2 ,fc},'**, and {tUA*,jt} are independent of the se- 
quences {w*} and {b k }. 

The direct way of evaluating £ and ?7 is to find the averages 
of e\ — a% and ||vjt|| 2 over a large number of independent runs 
for a fixed value of & in the steady state. An alternative simpler 
way, that can be used under the assumption that v k is ergodic, 
is to evaluate £ and 77 from one run by averaging over a large 

*ifimKa«> r\f ttAmt**\ne in tVia ctoorlir etot*» Qi-nr*** *i/a till VP s IflttTP 



TABLE II 

Analytical and Simulation Results tor 
A 7 = 2.cr b = 0.2, £ = 0.75. <ru- = 0 01 





CASE 1 


CASE 2 


Analytical 


Simulations 


Analytical 


Simulations 


-RLS 
Sain 
JLMS 
Snin 


0.5 


0.575 


1.325 


1.12 


Snin 
-LMS 
^rain 


1.25 


1.075 


1.25 


1.12 


RLS 
%in 

LMS 
\in 


0.8 


0. 86 


2 


1.77 


Vn 
LMS 
Xin 


1.2S 


1.15 


1.25 


1.18 



TABLE III 
Analytical and Simulation Results for 
N = 5,ff fe = 0-5, P - 0.75. «r H - = 0.005 





CASE 1 


CASE 2 


Analytical 


Simulations 


Analytical 


Simulations 


-RLS 
Siin 
pLMS 
Sin 


0.512 


0.68 


1.8 


2 


^in 
-LMS 

*min 


1.25 


1.43 


1.25 


1.2 


RLS 
Vn 

LMS 
\in 


0-73 


0.81 


3.13 


3.73 


SA 

LMS 
Vn 


1.25 


1.18 


1.25 


1.35 



Simulations have been carried out at several degrees of 
correlation between W2,k*, - • - . and w^k- Here, we give 
the results of the two following extreme cases: 

Case 1: uy+i,* = 3 = h 2, • ■ ■ , N - 1, for all k. 

Case 2: vy+i,* = Wj,k, j = 1, 2, • - * , AT - 1 , for all fc. 
Table II shows analytical and simulation results of a 2-D case 
(N = 2) while Table III shows the results of a 5-D case. In 
Table II, a h = Q.2,a w = 0.01, while in Table III, tr b = 0.5 
and <r w = 0.005. In both tables, p = 0.75. As seen by 
Tables II and III, both analytical and simulation results agree 
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' VIII. Conclusions-^ 

The tracking capabilities of the RLS, LMS, and sign algo- 
rithms for adaptive FIR fillers are compared in the case when 
the observation vector, x k , is stationary while the optimal 
filtering vector. r fc , is randomly time-varying. The following 
conclusions are drawn from the paper: 
(^))The steady-state excess mean-square estimation error, 
does not depend on the eigenvalue spread of the 
covariance matrix 7? = E(x h xl) in the cases of the 
LMS algorithm and the sign algorithm, while it does in 
^ the case of the RLS algorithm, 
r 2)\For the three algorithms, the steady-state mean-square 
deviation r/ between the weight vector of the adaptive 
filter and r k depends on the eigenvalue spread of JL 
(3) jFor the three algorithms, the value of the adaptation 
"^parameter that minimizes £ is different from the one 
that minimizes 7). 

4) The ratio of the minimum value of £ attained by the 
RLS algorith m to the one attained by the LMS algorithm 
is equal to tr (QR) /(tr (R) tr (Q)), with N being 
the number of adaptive filter weights and Q being 
the covariance matrix of the optimal vector increment 

5) The ratio of the minimum value of 77 attained by the RLS 
algorithm t o the one attained by the LM S algorithm is 
equal to y/tv(Q) tr (IT 1 ) f(Ntr (If 1 ))- 

6) A sufficient condition for the tracking of the RLS 
algorithm to be worse (resp. better) than that of the LMS 
algorithm is that the optimal vector increment c fc+1 — c fc 
be biased, in mean-square sense, toward the directions 
of the eigenvectors of R associated with the large (resp. 
small) eigenvalues. 

7) The tracking capabilities of the RLS algorithm are the 
same as the ones of the LMS algorithm in any one of 
the three following cases: (1) if R has equal eigenvalues, 
(2) if Q has equal eigenvalues, or (3) if R is diagonal 
and the diagonal elements of Q are equal. 

8) For Gaussian signals, the tracking performance of the 
sign algorithm is 1 dB lower than that of the LMS 
algorithm, independent of R and Q. 



^Appendix 

^EgOOF OF (6.4) AND 

With A,- and v t being as defined in Section VI, one has 



^rf- (U) 
Since tr(QR) = E^Rd^), then (LI) and (53) yield 

A* 

tr = £ VK- (1.2) 



Equations (I.I) and (6.3), respectively, imply that lr(R) = 
Efli Ai, and tr(Q) = YlLi 7i- Then, due to (1.2) one has 

A r N N 

AT tr (QR) - tr (Q) tr (R) = £ \ i7i - £ £ X ^ 

A" A" 

= E £ <&> 



where 



(1.4) 



One has 

A T A' A r A*-l A" 

E E *>v = E + E E ^ + < I5 > 



i-i 3=1 



Due to (L4) f Zij = Q.zij + z M = (A,- - A,-) (7* - 7>), and 
then (1.3) and (L5) yield (6.4). 
To prove (6.5), we use (IJ) and (6.3) to obtain 



N 

tr(g7?- 1 ) = ^A- 1 7i 



(L6) 



and then 



Ntr(QR- l )-tr(Q)tT(R- 1 ) 

N N N 

^E^-EEVV 0-7) 

Continuing in a procedure similar to that given in (I.3)-G-5), 
one can easily show that (1.7) implies (6.5). 
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Unitary ESPRIT: How to Obtain 
Increased Estimation Accuracy with 
a Reduced Computational Burden 
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Abstract — ESPRIT is a high-resolution signal parameter esti- 
mation technique based on the translation;*] invariance structure 
of a sensor array- Previous ESPRIT algorithms do not use the 
fact that the op erator representing thg^jiase_ddaysJietween the 
two aubarrays is unitary, fjfere . ^xejM^aentjij^^ 

cirdeTTfcenlr^sym^^ Unitary 
IcSPRlTrtB^^ ESPRIT- 
!ike structure except for the fact that it is formulated in terms of 
real-valued computations thr oughout Since the dimension of the 
matrices is not increased, this completely real-valued algorithm 
achieves a substantial reduction of the computational complexity. 
IBrthermore, Uni tary ESPRIT incorporates forward-backward 
a veraging , leading to an improved performance compared to 
$e standard ESPRIT algorithm, especially for correlated source^ 
Igials. Like standard ESPRIT, Unitary ESPRIT offers an in- 
expensive possibility to reconstruct the impinging wavefronts 
^jkggL£2E £>- These signal estimates are more accurate, since 
ffiBHry ESPRIT improves the underlying signal subspace esti- 
mates. Simulations confirm that, even for uncorrelated signals, 
the standard ESPRIT algorithm needs twice the number of 
snapshots to achieve a precision comparable to that of Unitary 
SPRIT. Thus. Unitary ESPRIT provides increased estimation 
accuracy with a reduced computational burden. 



Li Introduction 

THE recovery of signal parameters from noisy observa- 
tions is a fundamental problem in (real-time) array signal 
processing. Due to their simplicity and high-resolution capa- 
bility, ESPRTT-like subspace estimation schemes have been 
✓attracting considerable attention. Their parameter estimates are 
obtained by exploiting the rotational invariance structure of 
the signal subspace, induced by the translational invariance 
structure of the associated sensor array. This can be achieved 
without computation or search of any spectral measure [15], 
[17]. -Unitar y ESP RIT achieves, even more ,_aocin2fejPSSults 
ftan- previous ESPRIT^^ 

unit magnitude propert y of the phase factors that representjhe 
phase delays between the two suharrays [4]. It has been shown 
mTl2] that constraining the phase factors to the unit circle 
can also give some improvement for correlated sources. £ox 
centro-svmmetric sensor ar rays with a tran slational invariance. 
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structure, Unitary ESPRIT provides ^ very simple aad_e fficie nt 
solut ion t o this t ask. 

Alt hough Unitary ESPRIT effectively doubles the nu mber 
oj[j2fdata samp les, the comput ational complexit y is reduced 
by_j ransforming the requ iredLm nk-revealing factorizations o f 
complex mat rices into de j^npositio^ m atri- 
ces ot tho^ame sire. Thus, we obtain increased estimation 
accuracy with a reduced computational load. This reduction 
can be achieved by^cpjismjfiMgjn^^ tTaniformations that 
mar> centro-Herjcni.aan matrices to real matrices. These trans- 
formations have been introduced in Lee's pioneering work on 
centro-Hermitian matrices [10]. More than a decade later, her 
results were used to transform the complex covariance matrix 
of a uniform linear array (ULA) into a real matrix of the 
same size [8] to reduce the computational load of adaptive 
beamforming schemes [9], In this paper, we use more general 
centro-symmetric array configurations that have been receiving 
increased attention lately [22]. WL^lt^^^?^^ n L^^^ 
root version of Unitary ESPRIT that only requires reai-valued. 
computations from start to finish, by operating ..directly on 
^the^^insteajdo it to obtain sample covariance 

matrices.. It is well known that benefits result from smaller 
matrix. conditioning numbers [15]. With inrinite precision, both 
strategies would be the same, whether eigendecompositions 
-or singular value decompositions (SVD's) are used. Finite 
precision arithmetic, however, is employed in practical ap- 
plications. Therefore, numerical issues like round-off error 
and overflow are potential problems to be aware of when 
covariance matrices are estimated. I n addition to this square 
root approach, we also descri be an jdten^ 
covanance^proac jkjr^^ be more 

In the presence of additive noise, the computation of an 
optimal signal subspace estimate requires an SVD or an 
eigenvalue decomposition (EVD), which is computationally 
expensive, since 0(M Z ) operations are necessary to update 
the SVD or EVD if a new sample vector of dimension M 
arrives. Therefore, a number of alternative decompositions 
have been proposed to estimate the signal subspace in a 
computationally more efficient way. Examples include the 
rank-revealing QR decomposition [1], [2], the rank-revealing 
URV decomposition [11], [18], or a new Schur-type method 
for subspace estimation [3], [9]. These approximation tech- 
niques are computationally more efficient and well suited for 
a parallel (systolic) implementation, but they involve a certain 
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loss of accuracy that can be compensated by combining them 
with Unitary ESPRIT, yielding not only improved estimation 
accuracy, but also completely real-valued algorithms. In [5], it 
is shown how Unitary Schur ESPRIT dramatically improves 
the performance of the new Schur-type method, an adaptive 
subspace estimation scheme with a computational structure 
and complexity similar to that of a QR decomposition, except 
for the fact that plane and hyperbolic rotations are used.Jnthis, 
case, thej^jquked rank-decisia 
oTsigj^^s^ 

are^jguightfeiw^ The fully real-valued Unitary ESPRIT 
concept can also be extended to spatially smoothed forward- 
backward estimation schemes [6], [13], [14] and is applicable 
to many other subspace estimation techniques (see [20] for 
an excellent overview). The results are comparable to the 
advantages obtained by operating in beamspace [24] without 
the necessity of converting the data from element space to 
beamspace. 

This paper is organized as follows. It starts with a review 
of the definition and basic properties of centro-Hermitian 
^latrices. These properties will be used to derive the real- 
vSed implementation of Unitary ESPRIT. A brief review of 
^standard ESPRIT algorithm is given in chapter IE. It can 
beCseen as a generalization of the matrix pencil method [7]. 
Chapter IV introduces t he Unitary ESPRIT concept for centro- 
^yfhmetric array structures. IrTSe^ all 
d^3e~r^uTre^^ 

ifM^ecompodtions of real-valued matr ices of the same size 
yilffingTarcom A new reliability test 

which is a substantial improvement of current high-resolution 
aficly signal processing and spectral estimation techniques, 
is^f resented in Section IV-C. Further simplifications of the 
a||brithm are derived in Section IV-D, before a summary 
ojSlUmtary ESPRIT concludes the chapter (Section IV-E). 
E^ally, computer simulations compare the performance of 
Unitary ESPRIT with that of the well-known standard ESPRIT 
algorithm (Section V). 



^ilj centrq-Hernqtian Matrices 

^ First of all, let us introduce our notation and review the 
definition and the basic properties of centro-Hermitian ma- 
trices that have been derived by Lee [10]. Throughout this 
paper, column vectors and matrices are denoted by lower 
case and upper case boldfaced letters, respectively. U p is the 
p x p exchange matrix with ones on its antidiagonal and zeros 
elsewhere 

1 

~ ; 1 



.1 



€ R pxp . 



Since i7j is a s 



i.e. 



where the overbar denotes complex conjugation without trans- 
position. , 

Centro-Hermitian matrices of size (p xlp fonn, a; p • f/- 
dimensional linear space oyer R [10]. To" show how centro- 
Hermitian matrices can be mapped to matrices with real 
entries, Lee defines left JT-real matrices in the following 
fa shion. 

{Definition ^) 101: Matrices C7 Xf/ satisfying 



tetric permutation matrix, it is involutorial, 
= J 1 With this notation, we can define centro- 



are Jeft, H-rml- 
The unitary matrices 



n/2 



In 



jln 
-jUn 
0 

0 



~jHn 



(3) 



(4) 




M ^ TZ^MU, 



[G n p Gn q ]ec p 



X2q 



4~ 



* _ t 

Hermitian matrices in analogy to centro-symmetnc matnces. 

^iDenm^np A complex matrix M € C px? Js_ calleii 
centro-Hermitian if 



is centxo I JfcrjD2iJiari. Thus, the matrix 



for example, are left IT-real of even and odd order, re- 
spectively . More left 17-real matrTcesTcan Be oftairiedTby 
post-multiplying a left IZ-real matrix Q by an arbitm^M&aL ^/ 
m&lixJg) eyjggunatrjx we are 

in a position to state Lee's main result, which establishes an 
automorphism between centro-Hermitian and realjnatrices. 

Let T a and U, f denote arbitrary^ nonsin - 
Lgflicjes^ fsize p x p and (jX q, respectively. 
Then, the bijective rnappin 



ma ps the setofaU_ ^ x a centrqJ 3eiinitian^^ , 
thftjtffljg f all real matrices of the same size. 

Xbisjtheorem can, for instance, be used to ca}cjdatejhe7 

siQgular_y_aiue^4e£^ j 

^gg^^g^Let M be centro-Hermitian, and assume that 
the SVD of <p Q (M) = Q?MQ n € U pXq is given by 
&q(M) = U^EyVj, where the matrices Q p and Q q are 
unitary as well as left IT-real. Then, an S VDi of M is obtained | 



where jheJeft^d^right-singular. vectors of M are left /7-reaL 
Ewgf: The first part follows from the unitary nature of 
Q v and Q q , the second from the fact that the singular vectors 
of a real matrix are real. ■ 
For future reference, we, consider an efficient computation 
of a particular transformatipn T( ).\lt transforms an arbitrary 
complex matrix X7 6 C^jinto a real px2q matrix, denoted 
by~T(Gj. Notice that for every matrix <g, the matrix ^ 



!§^q([g n p Gn q \) ^ 
~Q p [g n p Gn q ]Q 2q -Sl^m 

J^Recall that the SVD of a complex matnx is .unique up to a umtarv diagonal 



-IT 



1234 
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is always real according to Theorem 1. Consider the case 
w here the left VT-rea l matrices Q v and Q 2 , ; are chosen ac- 
c ording to (3) or (4). Fu rthermore, let G be partitioned as 



G 



Go 



where theJiLodt^natrices & and Go should have th e same 
size. Qbj2<^&Jh^^ 

gjgn- Jnet^ ttgightforw aid calculationTshow that the desired 
real- valued magxx (6) can be ex p ressed as 



~9 




Re{Gi 4-ZTGo} 

V2*Re{/} 
Im{Gi + UG?} 



-Im{Gi - /TGo} 
ReJGx-ffGo} 



(7) 



[TANDABO-ESERIX 



^y^Stan dard ESPRIT Scenario . 

m Consider the standard ESPRIT scenario [15], [17], i.e., an 
i^M-element sensor array composed of m pairs of pairwise 
~ = identical, but displaced sensors (doublets). Let A denote the 
--d istance between the two subarrays. Incident on both subarrays 
C^are d narrow-band noncoherent planar wavefronts 

TO \s k (tn) = ujt^^^^l l£_fe<rf < m 

Qwith signal propagation velocity c and a common center 
^frequency wo- The d impinging signals are combined to a 
signal vector s(t n )~ Assume fo^^ 

subarrays donot share any elements, Le., they_do not overlap. 
T^enT the total number of sensors equals W =~2inj and the 
^ uncorrupted signals received at the two subarrays Jiaye the 
following form; ~ 



X(tn) = 



*l(*n)" 




A ' 


^2 (*n)_ 







S(tn) = AGs(t n ). (8) 



(Ag}€ £ M * d and6l£ C mX<i are the stemng^matrice^of the 
whole array configuration (global array steering matrix) and 
the first subarray, respectively. NTntk^tl^ tfe> frh columns 
of _ both array steering matrices depend on the .d irection of 
arrival (DO A) 8h ofjtheAt h source relative to the djsB laggnaent 

^be^em^ the_Lwo-subafrays. 2 Furthermore 

is a diagonal matrix of the phase delays between the sensor 
doublets for the d wavefronts. Its diagonal elements, the phase 
factors 0jfc T are given by 





Planar array co mp osed of m — 3 pairwisejdgtiticai but displaced 
seusorsj doublets ) . 



Here, Re{-} and Imf} denote the red and the imaginary 
part, respectively. Once again, if p is the center row 
of (7) should be dropped. Than* a n efficient computation of ? 
TiG) g W x2q from the complex matrix G^onlv, remiire^ 1 





Subarray 1 m 


= 5 




; ; i ; 

i I i . ; . 


? 




Subarray 2 






(a) 






Subarray 1 ,m 




1 


t r 
i ? 






Subarray 2 






(b) 






Subarray 1 m 


= 4 




t r M 


• 



Subarray 2 
fc) 

Fig. 2. Three different subarray choices for a uniform linear array (ULA) 
of \I = 6 identical sensors, (a) Maximum overlap (m ~ 5); (b) interleaved 
(m = 3); (c) mixed (m = 4). 

C& )More Structured Array Geometries 

Recall that every row of Ac corr esponds to an element 
of the sensor array. Tn the^ase_of (overlapping subarrays] 
a particular subarray con figurationi s described b y'seiection 
matrices that choosem elements^of x(t n ) e Q^Z^^^ 



m<L^^sJ^^ in each subarray . Let ^ 

and [J2J be (m x j?) selection matrices that assign elements 
of xffn) to the subarrays one and two, respectively. Fig. 2, 
for example, displays three different subarray choices for a 
uniform linear array (ULA) of M ~ 6 identical sensors. 

In generaL the^tw^ selectiOTLJB 
centro-s ymmetric with respegJLXo^one^jQiher^Le. 



a property that plays a key role in the derivation of the fully 
real implementation of Unitary ESPRIT, cf. Section IV-B. 
Therefore, the combined selection matrix 



Ji 

J2 



is centro-Hermitian, i 



Le.jffsri 



(9) 



© Qenerg ^ESPRJT Principle 
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are for med, obe ^ing^ 

JX=J[x{ti) x(t 2 ) 









A 1 












Abasis faote^stimated-agnaj-subspace-is-detennined from, 
the d dominanXlejBt.singul^..vectors.-according to , 



\JU„ =, 



Itiseasyjojeejfauev^^ 

qfjhe^ensjorjanry. l^ajionjn)jmjgto^^ 

X are rankndefident namely rank Xi = rank X 2 = rank 

T"^lTThus, ^e.cLcolumns_of- 



Pco! = 



A 



SP 



col 



form a b^js^ f^A^coJumn space of Xjf 

'SPaA € GL(d)7; (12) 

Here, GL(d) C C*** denotes the general linear group of all 
ioaiSijE) matrices of dimension ^Txj?5 Observe that^the. 
column space or range^of JX ^range^JX j: Ci"l,.is Ji&uaUy 
<^^^j^^!/^sggce A hi the same way, the d rows of 

[r, r 2 ]=p,ow[Ci c 2 } 

fUm a basis, for. the row space of [C± C 2 ] if 



Then, a unitary_basis for the row space of [Ci C 2 ] can 
also be obtained by computing its SVD (total least squares 
approach). However, it is 4ess.expensiveMq use jP row = C x < 
which corresponds to the standard least squares solution of 
the overdeternuned set of equations 



(16) 



-1; 



followed by an eigendecomposition of ^ - r x 
Mu ltiple Invaria nce Stature 

Unitary ESPRIT is applicable to centro-symmetric ( array) 
configurations. A sensor array is called centro-symmetric [22] 
if its element locations are^sy^ej^^ 
centroid and thq^cQmgle^^ 
arethe^ame. T heir globa JLagayistemn^ 
satisfies 7*ZXl 



[T GL{d)l (13) 

Merefore r _the rank-reducing numbers of the matrix pencil 

J [j$ ^ XTi= Prpw4l* - Md)SPcol. 

Jke the diagonal elements of # (f^e Japtprs) and can be 
Mfculated as the generalized . eigenyajues„ot^e_matrix pair 

IH^ac t( > *ese observations, the ESPRIT algorithm reduces 
^ Jd choosing the appropriate(^mpres^ that define 

Se required bases. In the absence^ noise (thecase discussed 
tsb far), any matrkes^^and^rW satisfying (12) and (13) 
will do the job. Withn^sy measurements, however, we are 
faced with the problem of estimating the signal subspace and 
'^s dimension. 

1 D.l SVD-Based Subspace Estimate 

The most robust way to estimate the required bases is to 
compute the singular value decomposition (SVD) of 




(17) 

f^sjomeoniitarxd^^ Notice that the 

matrig^j qA^ 2 is left JlreaL 

Uniform linear arrays, for e rn m pte t*"» m ™ f rnTr>TTton arrflVS 
used in practice, a re centro-symmetric. It is well known} 
that the analogy between array signal processing and time ^ 
series analysis (harmonic retrieval) can be obtained through ^ 
uniform linear arrays (ULA's) by interpreting them as uniformv 
sampling of a time series [16]. ^ 
Jhe ' centro-sjfflmetry of the global sensor array lA<i and 
imply Sat the steering matrices of both subarrays are 
a^^e^o-syjn«Retricr^ 

yeJMtarxJ^ matrix 

z d ={x n M x] 




x = p} 



u 0 



S 0 



r 5 H 



(14) 



admits the factorization 
JZ = 



where X denotes the measurement r^atrix X corrupted by 
additive, spatially uncorrelated 3 noise;\27) a®tmnsJts^Lfen^ 
inant singular values, and the unitary matrices U and V are 
partitioned accordingly. Then, the^esrirank£ar>g^^ 



X in the Frobenius-nonnTsgiven by ]X = U a £ 9 V? ) In other 
words, ^urBrrmr^timarofris the matrix X satisfying 

\\x-xy= l|Jr-y|| F . as) 

3\ rank Y<d 
*\ 

\ 

3 if the spatial covariance matrix of the additive noise is known up to a 
scalar factor, the SVD can beVeplaced by the generalized or quotient SVD 



X 2 



n m x 2 
n m Xi 



[S $~ l AS\ 



(18) 



whicfaj5-fiasily-seen by using-the .centro-symmetry of ihe 
subjirraxs_aiid_ the_uni©ry naturcjDf f. Thus, Z is also rank- 
deficient, namely 



show th ?T TTnitar y KSPR1T esse ntially doubles the nu mber o l 
available measurements from ALto^/V- Increased_«tin^tion 



accuracy^aruihejelore^eadi^ 

m^STSattSTx € C MxjV ^Q£the_standardESERrX formulation 



'-^-^^x^^igjjcojjespoiidsJ^ 
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f& Implementation^ 

Due to the special algebraic structure of the noise-corrupted 
data matrix Z and the structure of the subsequent total 
least squares (TLS) problem, the computational complexity of 
Unitary ESPRIT can be reduced significandy. This is ach ieved 
by jransforming the three f complex-value dljr ^c-reveali ng 
fact orizations, 

• t he subspace estimation step 

• thejubs equent total least squares problem 

• theJfoaL^en^^ 

into factorizations of real-valued matrices of the same size. 
Tnu5C"rea^ 

steps of the Unitary ESPRIT algorithm. The following three 
propositions derive the required transformations by taking 
advantage of the mapping between centro-Hermitian and real 
matrices, cf. Section H. In Remark 3, we also show how the 
real-valued total least squares problem can be replaced by a 
real-valued least squares (LS) problem. 
TS^q^^ Th e principal 

rjubspace of Z £ C Mx2iV (and, therefore, also the principal 
^ubspace of JZ) can be obtained through a rank-revealing 
^factorization of the real matrix T(X) G R Mx2N , where 
y[he transformation T(*) is defined in (6), Then, the complex 
Lmatrices Cx and C 2 , spanning the estimated signal subspace, 
%pbey 



Co — HmCil 



(19) 



\ frooQ By post-multiplying the noise-corrupted matrix 
\Z with a unitary permutation matrix we obtain a centro- 
5 Hermitian matrix in the following fashion: 



[x n M xn N ]- 



(20) 



According to Corollary 1, a rank-revealing factorization of 
Z C n can, thus, be obtained through an S YD of the real matrix 



(21) 



which proves the first part of this proposition. Let the d 
dominant left singular vectors of <pq{Zch) be denoted by 
E s € R llfx<£ . Then, the d dominant left singular vectors of 
Z C h as well as Z are given by Q M E S , Therefore, the matrix 



[£]-[: 



Jz 



provides a basis for the estimated signal subspace. With (10) 
and the left ZT-realness of Q M we, finally, get the desired 
result: 

— IT m i7 m J 2 /7 mQm&* 
= JzQmE* = 



Ero^osition l^Total Lej^ Sq^res Problem: The complex- 
valued SVD of size mx2d that solves the total least squares 
(TLS) problem C\& ^ C>, which is associated with Unitary 
ESPRIT, can be transformed into an SVD of the real matrix 
T{Ci) e R mx2t *, where the transformation T( ) is defined 
in (6). MoiSQff er, the eigeny Mie„s_j£fc_of the resulting TLS 
solution ^tls € C dXt£ will be symmetric with respect to the 
vuWlnrclen^ tHeSTafe Indices 



kJe{h 2 V - - : , dX^such that 



i 

f H S 



.(22^ 



^froafc^ The multidimensional TLS problem C$ « C> 
can be solved through an SVD of 

\2i 



[Ci a 2 \ = \u 1 r/,i 



I7-, 



Then, the TLS solution is obtained from the r/jight singu- 
lar vectors corresponding, to the d smallest singular values 
according to 



f(23 



*TLS = -^12^22 

where^we. have„ assumed that V 2 2 € GL(d). i.e., the TLS 
solution is unique. For the singular case, the reader is referred 

tM2ir~~~~ 

Thus, the TLS problem associated with Unitary ESPRIT 
can be solved through an SVD of 



[Ci 



n m Ci 



Notice that this matrix has the same structure as 

Z = [X 27.V/J]- 

Using, therefore, the same reasoning as in (20) and (21), the 
TLS problem is solved by computing an SVD of the real 
matrix 



r(c 1 ) = Q2[c 1 n m Cx] 



id 



Its right singular vectors will be denoted by 



W n W i2 
W 2 i w 22 



E U 



2dx2d 



Then, the right singular vectors of [Gi C 2 ] are determined 
from 



V s = 



v 12 v 22 



= W H Q? d 



(25} 



and ^tls is obtained from (23). Since the matrix Q^jW ls 
left iT-real, it can be written as 



for somejnamxZi-JE C dxU , cf. (2). With (25) we. therefore, 
conclude V22 —Yx^l Thus.jf <!>i is an eigenvalue of the TLS 
solution i^tls € GL(d),l/<£( is an eigenvalue of 
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The eigenvaU 

ues of the co3mpJex-jmttix^xLs _ c ^ determined from 
the eigenvalues" of a real matrix of .size d x d via the linear 
fractional transformation^ 



g + i 1 



(26) 



M oreover, th e eigenvectors of both matrices are Jdentical. 

^\^£o^(^} Assume, for the moment, that the left JT-real 
matrix Q 2d is the one we have defined in (3). Then, (25) yields 



In 
In 



Jin 

~jln 



w. 



After partitioning V and W as before, we therefore conclude 
from (23) 

*TLS ="(^12 +jW 22 )(W l2 -jW 22 r l 

^ ^-((-WnW^-jlja-WnWg) + jh)' 1 
i -/(^tls) with T TLS = -W 12 W 22 K (27) 

Sfere, f{x) denotes the linear fractional transformation (26), 
^Bich is analytic for x ^ — j. Let 



Proposition 2 states thatJthe eigenvalues of #tls? i-e^ ^ i 
phase facto rs estimatedj ja Unitary ESPRIT, are symmetric 
with respect to the unit circle, since they satisfy (22). This 
observation gives rise to a new reliability test provided by 
Unitary ESPRIT without demanding additional computations. 
This reliability test is a substantial improvement of current 
high-resolution array signal processing and spectral estimation 
techniques since usually there is no easy way to determine^ 
how reliable the resulting estimates are.^Unr^Ii^jjle^^ 

^si^tsmight have been cau^d^hy^jfaise_^e&tjipate„o^hej 
number of sources ^f"or""6y theJgcLik^tbefe^s np^ source 

f 

^^Eitnar k 2p ■ Ei genvalues with Unit Modulu s^ Notice first 
that the eigenvalues <j)k that lie on the unit circle form a subset 
with nonzero measure in the class of all eigenvalues fulfilling 
(22), i.e., being symmetric with respect to the unit circle. 
Owing to this and the fact that Unitary ESPRIT produces 
consisten t DO A estimates, asymptotically all ihe , estimated 
phasej actors 6u will be on the unit circle. 

If, however, the number of snapshots N is too small or if 
there is only noise present, the eigenvalues of ^tls might 
fail to satisfy 



Ttls = ~W l2 W£ = TOT' 1 



(28) 



w Q 2 d=<?2d^2d where R 2d E R 2dx2d . 

After replacing the real matrix W by W = R 2 dW. we invoke 
the same reasoning as above to prove this proposition for an 
arbitraryj eft JT-real t nmsformatio n6 2 rf • ■ 
j ^gmgr&il -^^ Ap pmachl Instead of the described 

square root (or direct data) approach based on a real-valued 
SVD, cf. Proposition 1, we can use a covariance approach 
based on a real-valued EVD to determine the signal sub- 
space estimate. Then, B 9 £ R Mxd denotes the d principal 
eigenvectors of 



T(X)T(X) H e U MxM . 



(29) 



First, forming T (X) according to (7), followed by the compu- 
tation of (29), is more efficient than the approach alternative 
suggested in [8] and [24]. There, it is proposed to compute the 
complex- valued sample covariance matrix Rxx = XX € 
qMxm first Th cn ^ i s determined from the EVD of 
RelO?,flvvOwV which is comnutationallv mnrp axnMKive 



|<M = i 



(30) 



bg the eigenvalue decomposition (EVD) of the real matrix 
Ttls * It is a well-known result J^m^fuaction^theory- that~the- 
<gg envaLues-ofjS^Lsjcan be o btained through the samejinear 

fractional transformation, i.e. ~ 

— , — „ — — — ^ — 

\ y $ = witft & = diag{a?fc}fc =1 and u k ^ -j 

agd the conesponding^i^nv^tors^oXXxLS^^d ^TLSjare, 
jL4pntical. ~ 

vb) ;An arbitrary left XT-real matrix of dimension 2d x 2d 
can, obviously, be written as 



which indicates that the subsequent estimates will be unreli- 
able. Hence, no further computations should be carried out. 
Condition (30) implies that all eigenvalues Uk of T T ls are 
real, cf. (26). 4 jTht^ , ijLsomej>fjte 
^ojuLgatejgirs, theUnitary& 

^d L giealionthm h^ S^£55arted^^ an increased window 
length jV" or more reliable measure ments. If, conversely, all / 
eigen^l^s^ Aare real, i.e., the reliability test has been \ 
passed," alTbsnmated phase jactors <f>k are precisely on thev 
it circle. 



P.) Real-Valued Least Square s 

Notice that the derivation of the Unitary ESPRIT reliability 
test is based on a total least squares solution of (16). Thus, 
the computation of Ttls requires an SVD (or another rank 
revealing factorization) of T{C\) £ U mx2d . By computing the 
less expensive, least, squares instead of the total least squares 
solution of (16), we would, however, lose the benefits of 
the reUability test, since j(22X would , no, jonger be. satisfied. 
Moreover, the complexwalued least squares problem (16) 
cannot be transformed into a real-valued problem of the same 
size, T^gEawi^ 

valued sauaresjwb^m^^blcli ^aaJ^^olyecLinstead. < 

^emar^3^Least Squares Estimate : After partitioning the 
reanSatrix of (24) according to 



nCJ^lTx T 2 ] 



with T U T 2 E 
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it is easyjo s ee tha t Ttls is a TLS solution of the real- valued 
system of equations 

To save computations, we can, therefore, solve (31) by com- 
puting its least squares solution Tls- Here, the Unitary ES- 
PRIT reliability test is still applicable, since the resulting 
matrix Tls is always real. If the reliability test has been 
"passed", the estimated phase factors are on the unit circle. 

The real- valued LS or TLS problem (31) can directly be 
obtained from E H by observing 



Now, we are in the position to summarize the described real 
implementation of Unitary ESPRIT, which is given in Table 
I. Here, th e left 27-real matrices Q m and Q M are chosen 
acco rding to (3) or (4). 

atice that a li near estimat e of Jhe source signal matrix ^ 
^StegjQ-£!^^ 

thi^^tionjojtf^^ 

J^J^^hejoeJ^ithoitf^^ 

^umedjoj^e-spatial^^ 



. Computer Simulations 



t(c 1 )=qZ[c 1 cyp Ed \Q 2d i£ 

= —^[KiEs K<>E<\ 

where the selection matrices Ki and K 2 are defined ja$ 
f ollows:^ ~~ ~~ 

•« \ k 2 =Q%j(Ji - n m j x n M ) QtJ 

= \. - i _ 

= Since th e matrices in braces are centro-Hermitian, Ki and 
I Ko are alway^jeaL cf. Theorem L They are even sparse if 
Jhe selection matrix Ji is sparse and the matrices Q m and 
iJQ M are chosen according to (3) or (4). This is illustrated by 
"the following example. For the-ULA ^ overlap 
^^tehedjin(^g. 2(a)) J^i^giwnJ^ 

?1 ~~ rl 0 0 0 0 0 

0 1 0 0 0 0 

0 0 1 0 0 0 

0 0 0 1 0 0 




<32)j 



Lo o o o i oj 



Thus, stnught forwanT calc ulations yield 

i i o\ [o~T~o|- 
a „ 1. 1 \ \CL.o_oi 

0 0_ js/2 0 0 0 

etT oT IP; i l o^ 
& gl.JP> oi l. 
ro o 



0 
0 



0 
0 
0 

-1 



-1 
6" 

0 
0 



-1 

T 
o 
o 



0\ 

ij 

0 

0 J 



In this section, we present some simulation results that com- 
pare Unitary ESPRIT with the standard ESPRIT algorithm, 
using the SVD implementation in all cases. Among others, we 
examine scenarios where the standard ESPRIT algorithm faces 
some problems, like l ow signal-to-nojs ejauos^ shojtwindow 
lengths, and con^te^qu^ 

s^jSignal Reconstruction 

First, we examine the effect of Unitary ESPRIT on the 
resulting signal estimates. To this end, three impinging wave- 
fronts are reconstructed using a single ULA of M = 9 sensors 
with maximum overlap, cf. Fig. 2(a). The three uncorrected 
equi-poweretLOP SK signals arrive from 9 X = JL0° f 9% =20? , 
andj?3,=_3Qf\ resr^ctiveiy. Fig. 3 depicts the resulting output 
signal-to-noise-and-interference ratio (SNIR) as a function of 
the SNR and the number of snapshots N using standard 
ESPRIT (dashed lines) and Unitary ESPRIT (solid lines). The 
values of N marked on the right side of the figure correspond 
to the solid lines, i.e., Unitary ESPRIT. The output SNIR 
achieved by the standard ESPRIT algorithm for a given value 
of N (dashed lines) can be found below the corresponding 
solid lines. For small values of N, e.g., N - 5 snapshots, 
U Unitary ESPRIT achieves a significantly better performance 
than the standard ESPRIT algorithm. Notice that standard 
^ESPRIT with N = 10 snapshots attains the same performance 
3 .as Unitary ESPRIT with N = 5 snapshots for SNR's that 
are greater than 15 dB, while the performance of standard 
ESPRIT with N - 20 is comparable to the performance of 
Unitary ESPRIT with N = 10 for SNR's that are greater than 
5 dB. TTju^Unit^ doubies,the~nurnber 
ofj^ilablejtnapshots JSL compared to the standard ESPRIT 
algorithm. 




E.j^ ummary of the Algorithm\ 
Beforejggsentjr^^ ESPRIT , we jiote. 



B. )QO A Estimation , 

Next, we investigate the effect of Unitary ESPRIT on the 
estimated phase factors ^ 1 < A; < d. Consider a ULA with 
M = 6 sensors and three correlated signals .impinging from 

^~~and <9 3 == 20° 



an interesting relationship brtwen ftigff^v^hi^ tf™ ™* 
rn ^x^<fe R Qted by and the^esrimated phase factors 
cf. (9). Solving <gJ>* = /(wfc) for the spatial 

fregueoeks^ 



is given by 



Their correlation matrix 



1 p p 

PIP 
P P 1 



(34) 




2arctan a>jfc, 



k = 1,2, 



The phase factors 4>u*hi and <£ 3 , estimated with the standard 
ESPRIT algorithm and Unitary ESPRIT, are marked by crosses 
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r JG&BL&uD 

f Summary of t Intta gv ESPRIT 



Q}. Initialization: Form the matrix X 6 C u x N from the available measurements. 

(2\signaISubspace Estimation: Determine tnc real matrix T[X) € R 6 ** 2 * ta^fgjand compute the SVD 
of T(X) (square root approach) or the cigendecoraposition olT(X)T{X) H (co variance approach). The <i 
dominant left singular vectors or eigenvectors will be called 6 R** - . Estimate the number of sources d. 
if d is not known a priori [22]. 

Ci^jTotal) Least Squares: Solve the overdetermined system of equations 



*r— — - - — : 

by means of least squares (or rota/ least squares) techniques. The selection matrices K\ and K2 are defined 
\4) Eigenvalue Decomposition: Compute the eigendecomposition of the resulting solution 



Y = T fl T" 1 € B** - , where fl - diag{w fc }J =1 . 



^Reliability Test: Ifail eigenvalucsurk are real, the estimates will be reliable. Otherwise, start again with more 
measureinenis. 

^OOA Estimation: Fsrimarr the directions of arrival (DOA's) from 



1 



/** = 2arctanid*| fc = 1,2,.. .,d^ 



according to (9). L ^ ^ 



where D f C Jx<i de notes an arbitrary diag onal frow> scaling matrix 141- 



»1? 




0 5 10 15 20 25 30 

SNRindB 



Fig. 3. Output SNIR as a function of the SNR and the number of snapshots 
N using standard ESPRIT (dashed lines) and Unitary ESPRIT (solid lines) 
for $1 = 10°, 0 2 = 20° T and 0 3 = 30°{.V = 9 sensors, 1000 trial runs). 
The values of N marked on the right side of the figure correspond to the 
solid lines, i.e., Unitary ESPRIT, The output SNIR achieved bv the standard 

FSPRTJ ^igorithm-far i ipvrn initio nf V f rin nhnH Ti nas) r^ r j h ft fpu m? hrlayr 
the corresponding solid lines . 

a correlation coefficient of p = 0.5. The results of 80 trial 
runs with N — 20 snapshots and an SNR of 0 dB are 
shown. Notice that all phase factors estimated with Unitary 



Standard ESPRIT 




Fig. 4. Phase factors d>i,&>r and 03. estimated with the standard ESPRIT 
algorithm for &\ = -20°, 82 = 0°,#3 = 20°. and correlation coefficients 
Pv2 = 0.5. P13 = 0.5, and m = 0.25 (M = 6 sensors, SNR - 0 dB, 
:V = 20, SO trial runs). 



coefficient of p = 0.8. In this example, the Unitary ESPRIT 
reliability test has failed three times. To picture these failures, 
the corresponding phase factor estimates are surrounded by 
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Unitary ESPRIT 




270 



wrig. 5. Phase factors oi*<h.i and 03, estimated with Unitary ESPRIT 
LJfor $\ = -20°, 0* - Q°.6>3 = 20°. and coiTeiation coefficients 
KUP\2 — 0.5. pi 3 = O.o, and p?3 = 0.25 (M = 6 sensors, SNR = 0 
^ dB, y = 20, 80 trial runs). 



Standard ESPRIT 




fig. 6. Phase factors 01, 02, and 03, estimated with the standard ESPRIT 
algorithnufor $1 = -20°, # 2 = 0° , 6^3 = 20° , and correlation coefficients 
p l2 =<0-8,)>t3 =(0.S,\nd P23 =(0.64W - 6 sensors, SNR = G dB. 
*V = 207&) trial runs£ " v — ^ 



lower than the variance of the DOA estimates obtained by 
the standard ESPRIT approach. The advantages of Unitary 
ESPRIT become even more evident if the root mean squared 
error (RMSE) of the estimated directions of arrival is plotted 
as a function of the correlation coefficient p. Fig. 8 show 
these curves for SNR's of -3, 0, and 5 dB using 3000 trial 
runs. Hie standard ESPRIT algorithm (dashed line ) is 
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Unrtary ESPRIT 




270 



Fig. 7, Phase factors 01.02* and 03, estimated with Unitary ESPRIT 
for $1 — -20°, 9n =z 0°,#3 =s 20°. and correlation coefficients 
P12 = O.S.P13 = 0.8. and P23 ~ 0.64 (M = 6 sensors, SNR = 0 
dB, JV = 20, 80 trial runs). Estimates that produced a failure of the reliability 
test are surrounded by a circle (o). 




correlation coefficient 



Fig. 8. Root mean squared error (RMSE) in degrees of the estimated 
directions of arrival as a function of the correlation coefficient p and the 
SNR for 0i = -20°. On = Q<\ and 0 3 = 20°(M = 6 sensors, N = 20, 
3000 trial runs). The signal correlation matrix is given by (34). 

line - - *) ? and Unitary ESPRIT with the new reliability test 
(solid line — ). It can be seen that Unitary ESPRIT improves 
the estimation accuracy considerably. In the case of low 
SNR's, the estimation accuracy is improved even further, by 
exploiting the information provided by the new reliability 
test. The corresponding failure rates of the Unitary ESPRIT 
reliability test are plotted in Fig. 9. 

-D ue to the forward-backward averaging effect. U nitary 
ESPRIT can separate two completely c^ hejEmt-&ax^^ 
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0.1 



0.2 



0.3 0.4 0.5 
correlation coefficient 



0.6 



0J 



0,8 



0.9 



Fig. 9. Failures of the Unitary ESPRIT reliability test as a function of the 
correlation coefficients p for #i = -20°, # 2 = 0°, and H - 20° {M = 6 
msors, N = 20. 3000 trial runs). Once again, the signal correlation matrix 
V% given by (34). These curves correspond to the solid lines in Fig. 8. 



1.5. 



1.3- 



1.2 



0.9 



Standard ESPRIT (LS) 

Standard ESPRIT (TLS) 

Unitary ESPRrr (LS) 

Unitary ESPRIT (TLS) 
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0.2 



0.3 



0.4 0.5 0.6 
correlation coefficient 



0.7 



0.8 . 0.9 




200 

degrees 



^Exg^JJ^ RMSE (in ^greesJ_of_Uie_ estimated directions of arrival as a 
function of trie' magnitude and phase of the complex correlation coefficient p 
f6r#r=~8°~and ^0°JisingIstandard ESPRITi-W* = 4 sensors, X = 
207"lT)0"trial runs). " - 



Fig. 10. RMSE (in degrees) of the estimated directions of arrival as a 
function of the correlation coefficient p for 0i = 0° and #2 = 20°{il/ — 4 
sensors, N = 20, 100 trial runs). Notice that the curves for the LS and the 
TLS version of Unitary ESPRIT fall on top of one another. 

signals with correlation coefficient p are impinging an a ULA 
of M = 4 sensors from 9 1 = 0°and 9 2 = 20°. Fig. 10 
shows the resulting RMS error of the estimated DOA's as 
a function of p. The performance of Unitary ESPRIT is not 
effected by the correlation, while the performance of standard 
ESPRIT deteriorates dramatically as p increases. Notice also 
that the difference between TLS and LS version of standard 
ESPRIT is negligible, while the LS and the TLS version 
of Unitary ESPRIT fall on top of one another. JThus, it is 
advisable to use_ t he LS version of Unitary ESPRIT instead 
of the computationally more expensive TLS version . Finally, 
Figs. 1 1 and 12 show the RMS error of the estimated DOA's 
as a function of the magnitude and phase of a complex- valued 

rnrri»latirni ropffirif nt n rnnfirmirKT tht*- rnnHncinnQ drawn 




200 

in degrees 



Vgg. ,12; \RMSE (in degrees). of the estimated directions of arrival as a 
function of the magnitude and phase of the complex correlation coefficient p 
for #i = 0° and On = 20° using Unitary ESPRIT (M = 4 sensors, N = 
20, 100 tnal runs). * " - J ' " 



VI. Concluding Remarks 

An improved version of the ESPRIT algorithm, called 
Unitary ESPRIT, has been presented in this paper. Uni- 
tary ESPRIT represents a simple method to constrain the 
estimated phase factors to the unit circle, yielding more 
accurate signal subspace estimates. The computational com- 
plexity is reduced significantly by exploiting the one-to-one 
correspondence between centro-Hermitian and real matrices, 
allowing a transformation to real matrices, which can be 
maintained for all steps of the algorithm. Jjnitary ESPRIT 
also provides a new reliability test, which is particularly 



useful in extremely low SNKls . Duelo the inherent forward- 
backward averaging effect, Unitary ESPRIT can separate two 
completely coherent sources and provides improved estimates 
for correlated signals. Moreover, Unitary ESPRIT offers a 
great potential to improve the performance of approximate 
signal subspace estimation techniques, which are well suited 
for an adaptive implementation, since inexpensive updating 
strategies are known [5]. 

Thp fptrt rhnt FTnitflrv F^PRTT i« pfftHpTiriv fnrmiitafpH in 
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critically important for the extension to 2-D centro-symmetric 
arrays with a dual invariance structure. 2-D Unitary ES- 
PRIT [23] provides automatically paired source azimuth and 
elevation angle estimates along with an efficient way to re- 
construct the impinging wavefronts. Furthermore, an efficient 
DFT beamspace implementation of Unitary ESPRIT has also 
been derived in [23], enabling reduced dimension processing 
in beamspace, if there is a priori information on the general 
angular location of the DOA's. 
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Abstract — In code division multiple access (CDMA) mobile 
radio systems, both iittersymbol interference and multiple ac- 
cess interference arise which can be combated by using either 
elaborate optimum or favorable suboptimum joint detection (JO) 
techniques. Furthermore, the time variation of the radio chan- 
nels leads to degradations of the receiver performance. These 
degradations can be reduced by applying diversity techniques. 
Using coherent receiver antenna diversity (CRAD) is especially 
attractive because only the signal processing at the receiver must 
be modified. In the present paper, the application of CRAD to 
ffle more critical uplink of CDMA mobile radio systems with 
suboptimum JD techniques is investigated for maximal-ratio 
combining. The authors study six different suboptimum JD tech- 
ttiijues based on decorretating matched filtering, Gauss-Markov 
Estimation, and minimum mean square error estimation with and 
without decision feedback. These six suboptimum JD techniques 
^?hich are well-known for single antenna receivers are extended 
Mr the application to CRAD. A main concern of the paper is 
the determining of the SNR performance of the presented JD 
techniques for CRAD and the achievable average uncoded bit 
error probabilities for the transmission over rural area, typical 
arban and bad urban mobile radio channels are determined. 



fy I. Introduction 

TWO MAJOR PROBLEMS must be solved in digital 
mobile radio systems, namely the multiple access (MA) 
problem, arising from the simultaneous transmission of signals 
associated with several active users that share the same trans- 
mission medium, and the equalization problem, arising from 
time- varying and frequency-selective mobile radio channels 
[1], [2]. An attractive solution to the MA problem is code 
division multiple access (CDMA) [1] which is currently under 
discussion for the application to third generation digital mobile 
radio systems [3]-[7]. In a CDMA mobile radio system, 
multiple access interference (MAI) between data symbols of 
different users occurs which can be combated successfully by 
applying either joint detection (JD) [4], [8], [9] or interference 
cancellation (IC) [10] techniques. Furthermore, intersymbol 
interference (ISI) arises between consecutive data symbols 
associated with a single user when CDMA is used due to the 
frequency selectivity of mobile radio channels. ISI is negligible 
when the symbol period T s is much greater than the duration 
of the channel impulse response. However, ISI becomes severe 
when T s is lower than or approximately equal to die duration 
of the channel impulse response. With respect to generality, 
the existence of ISI in conjunction with MAI in CDMA mobile 

Manuscript received September 28, 1993; revised December 13, 1993. 
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radio systems is assumed in the following. Therefore, a major 
requirement for detection in such CDMA mobile radio systems 
is the combatting of both MAI as well as ISI, respectively. 

In this paper, JD techniques fulfilling the aforementioned 
requirement such as presented in [4], [8). [9] shall be con- 
sidered. Unfortunately, the maximum likelihood JD technique 
of [8} is prohibitively expensive. Therefore, suboptimum JD 
techniques [3], [4], [9}, [1 1]— [14] are more favorable. These JD 
techniques are appropriately modified versions of wet! known 
adaptive equalizers such as reviewed by Qureshi in his famous 
paper [15]. 

Many mobile radio systems such as the successful pan- 
European Global System for Mobile Communications (GSM) 
[16], [17], the Digital Cellular System (DCS) 1800 [18], the 
Japanese Digital Cellular (JDC) [17] and the American Digital 
Cellular (ADC), also termed IS-54, apply discontinuous, i.e., 
burst, transmission. These systems, therefore use time division 
multiple access (TDMA) combined with frequency division 
multiple access (FDMA). Burst transmission is advantageous 
with respect to the implementation of the receiver hardware, 
cf., i.e. [16], [17] for a detailed rationale. Therefore, it stands to 
reason to develop novel mobile radio systems in an evolution- 
ary manner by introducing an additional CDMA component to 
the aforementioned F/TDMA based mobile radio systems in 
order to achieve a capacity improvement over the aforemen- 
tioned F/TDMA based mobile radio systems. Such a CDMA 
mobile radio system using FDMA and TDMA and applying 
JD techniques which was evolved from GSM was proposed 
in [19]. This system shall be termed JD-CDMA mobile radio 
system in what follows. Following the report of [19], burst 
transmission shall be considered in what follows. With respect 
to a moderate hardware expense of the receivers for the 
JD-CDMA mobile radio system, synchronization between the 
users, i,e., mobiles and base stations, is favorable and can 
be easily achieved based on the TDMA component, i.e. by 
following the procedures used in GSM, cf. [20]. However, the 
problem of synchronization shall not be covered in this paper. 
Furthermore, an isolated cell of the cellular JD-CDMA mobile 
radio system environment is considered in what follows. 
Intercell interference is treated as additional noise. 

In the synchronous JD-CDMA mobile radio system, the 
signal processing for the downlink and for the uplink differ 
considerably. In the case of the downlink, the signals asso- 
ciated with K simultaneously active users are radiated from 
the same location, namely the base station. Hence, all K user 
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which are radiated from A' separate mobiles, is associated with 
A' separate radio channels. Therefore, the separation of the A' 
user signals is more critical in the uplink than in the downlink. 
The authors consider only the uplink in the following. 

Conventionally, matched filtering, in particular by applying 
decorrelating matched filters (DMFs), or correlation including 
RAKE concepts are applied in CDMA mobile radio systems 
in order to accomplish the suboptimum data detection [I], 
[7]. [21 J. Matched filtering is often inefficient because a tight 
power control is required and the MAI is treated as noise. 
However, the data detection with DMF's can be regarded 
as a JD technique which shall be done in what follows. 
Setting out from the DMF, decision feedback equalizers which 
are similar to interference cancellation techniques, cf. [10], 
[22], can be implemented. Such equalizers shall be termed 
decorrelating matched filter block decision feedback equalizers 
(DMF-BDFE) in what follows. Further well known subopti- 
mum JD techniques for single antenna receivers applicable 
to a synchronous JD-CDMA mobile radio system are the 
zero forcing block linear equalizer (ZF-BLE) [3], [4], [11 J, 
[12) based on the Gauss-Markov estimation [23] and the 
minimum mean square error [23], [241 block linear equalizer 
(MMSE-BLE) [3], [12]. Setting out from the ZF-BLE and the 
MMSE-BLE. block decision feedback equalizers referred to 
as zero forcing block decision feedback equalizer (ZF-BDFE) 
and minimum mean square error block decision feedback 
equalizer (MMSE-BDFE), respectively, were developed for 
single antenna receivers [3 J, [UJ, [12]. 

The time variation of the radio channels leads to a con- 
siderable variation of the signal-to-noise ratio (SNR) at the 
detector input. With respect to a constant expectation of the 
SNR, the error probability increases with increasing variation 
of the SNR [2], Hence, the reduction of this variation of the 
SNR, i.e by using diversity techniques [25], is a desirable 
asset. The application of receiver antenna diversity [25] is 
especially favorable because only the signal processing at 
the receiver must be modified whereas the signaling scheme 
of a synchronous JD-CDMA mobile radio system remains 
unaffected. In the course of the paper, the authors introduce 
novel representations of the JD techniques DMF t DMF-BDFE, 
ZF-BLE, ZF-BDFE, MMSE-BLE. and MMSE-BDFE for 
coherent receiver antenna diversity (CRAD) and maximal-ratio 
combining [25]. Furthermore, it shall be shown that all six 
JD techniques are special cases of a basic concept for JD. 
A main concern of the paper is the determining of the SNR 
performance of these six JD techniques for CRAD and the 
achievable uncoded bit error probabilities for the transmission 
over rural area, typical urban, and bad urban mobile radio 
channels by using the channel models defined by COST 207 
[26], [27] are determined. 

The paper is organized as follows. In Section II, the uplink 
model of the JD-CDMA mobile radio system with CRAD 
which is to be investigated is introduced. The extension of 
the DMF, the DMF-BDFE, the ZF-BLE, the ZF-BDFE, 
the MMSE-BLE, and the MMSE-BDFE from the single 
antenna case to CRAD as well as the corresponding SNR 



transmission over rural area, typical urban, and bad urbk 
mobile radio channels are given. \ 
In the present paper, vectors and matrices are in bold face. 
Furthermore, the symbols (■)*. (-) r , Jj • ||, and £{*} designate * 
the complex conjugation, the transposition, the vector norm, 
and the expectation, respectively. 

IL System Model 

The discrete-time lowpass representation of the uplink in 
a synchronous JD-CDMA mobile radio system with CRAD 
shall be described in this section. The basic structure of the 
uplink is depicted in Fig. L It is assumed that K users, i.e. 
the mobiles, are simultaneously active in the same frequency 
band transmitting the finite data symbol sequences 

d<*> - {dfK 4 k) • • • <#>) r , € V. V C C. 

of N m-ary complex data symbols cfn ] with period which ^ 
are taken from the complex set * ^ 

e£.^= i.---.m. rn€N. (2) 

The actual data rate may be varied due to the service provided 
or due to the desired transmission quality. Each of the data 
symbols dn \ n = 1. - • - .iV, of mobile k is spread by using 
the user specific signature sequence 

4 = 1.-. -.JSC. g=L- K.QeN (3) 

at the transmitters in order to allow a coexistence of A' 
simultaneously transmitted signals in the same frequency band. 
The m-ary complex signature elements c q k ^ of (3) which are 
taken from the complex set 

V c = {tV.i- tf c o * ■ • fern}- v c 4 c 6 C 

= 1.* - ■ .m, m e N (4) 
are termed chips. The chip duration is given by 



Each mobile isjisjumedj 

fhT^Smitted signajs are. received, M the uplink receiver^ 
i.e. the base station, pyeiJf^ receiver antennas. Therefore, the 
transmission of the K user signals takes place over K • K a 
different radio channelsLwith the time- variant complex impulse 
responses 

where the radio channel with the impulse response f^ kka ] {rJ) 
fefers^thTc6¥hection roFrhooTFeTTwith receiver antenna k a . 
In (6), r denotes "the delay parameter referring to the time 
spreading, i.e. distortion, of the transmitted signals due to 
multioath reception, and f denotes the real time referring to 
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Fig. I. Structure of the uplink of a JD-CDMA mobile radio system with CRAD. 



K • K a radio channels are assumed to be wide sense stationary 
uncorrelated scattering (WSSUS) multipath channels [26], and 
therefore the simulation models of these radio channels can be 
based on the COST 207 specifications [26], [27]. 
in the simulation programs, the K - K a WSSUS multipath 
|ehannels are implemented in accordance with the Monte 
^Carlo-based model introduced by Schulze [28] and further 
invejggated by Hoher [26]. This simple statistical model is 
entirely determined by the two-dimensional scattering function 
$(tSd) which is only dependent on the delay parameter r 
andyhe Doppler frequency f D associated with the motion 
of the. mobiles. In the following, a brief overview of this 
statfsjjcal model shall be given. In the statistical model used 
in tte simulations, the distribution of the Doppler frequencies 
fo jjg assumed to be independent of the distribution of the 
dela^ parameter r. Therefore, the scattering function can be 
displayed as the product of the Doppler power spectrum 
^(tjfc) which only depends on the Doppler frequency f D 
and ; jE§e delay power spectrum qt(t.O) which only depends 
on tJil delay parameter r, i.e., 



S(T./ D ) = S f (O./ D )-0r(T.O) (7) 

N 

holds. The K ■ K a channel impulse responses hS k - k ^{r. t) of 
(6) are samples of K - K a bandlimited stochastic Gaussian 
processes which are assumed to be ergodic. Furthermore, 
only Rayleigh fading is considered throughout the paper. 
In the digital simulations, each hS k ' k ^{rJ) is the linear 
superposition of P, P 6 N, uncorrelated echoes, each with 
a particular delay parameter Tp k ' ka \ p = 1.---.P, which 
are represented by P sinusoids with equally distributed phase 

angles 9 { p k ' ka) £ [0. 27r), p — 1 P 7 and Doppler frequencies 

Axp * * P = 1- - • * * P* being distributed according to the 
desired Doppler power spectra associated with the aforemen- 
tioned bandlimited stochastic Gaussian processes. The Doppler 
power spectra are assumed to be the classical Doppler power 
spectra [29, ch. 1.2.1] in this work. The delay parameters 
-p " t) are exponentially distributed according to the delay 
power spectra associated with the desired environment types 
rural area (RA), bad urban (BU) or typical urban (TU) [27]. 
With (*?(*) designating Dirac's delta-function, the K - K a 



channel impulse responses are therefore given by 
) (r. t.) = Urn -j= £ exp{ > } 

• exp{j2nf { ^k} .*(t- t^) (8) 

[26]. A phase modulation fading simulator can be realized 
by implementing (8), cf. [29, Fig. 1.7- 1 (a)]. The factor l/s/T 
ensures the convergence of (8) in the limit P — * x, [26]. 
The Gaussian distribution of the quadrature components of 
h {kMa) {r,t) defined by (8) follows from the central limit the- 
orem provided that P is sufficiently large. In the simulations. 
P equal to 600 is chosen. In (8), the amplitudes of the various 
sinusoids are all equal to 1 / \/P which is reasonable because 
the simulation model of the radio channels is statistically fully 
determined by the two-dimensional scattering function [26]. 

In Fig. 2, a typical simulated plot of the absolute value 
\h^ k - ka \r. t) \ for arbitrary k and k a is depicted versus r and t 
in the case of BU and a vehicle speed equal to 50 km/h. The 
delay parameter r is resolved with 0.30 /xs in Fig. 2. In order 
to relate the duration of h^ k ' k<t) (r. t) to the symbol period i;, 
r is given in units of T, with T s equal to 7 /xs which is used 
in the simulations discussed in Section IV. A new sample of 
\h {k - k ^{r.i)\ is sketched every 100 along the /-axis. The 
time t is depicted in the range 0 to 8 ms in ail cases. The time 
variation of M****)(t. t) is obvious from Fig. 2. The duration 
of /^*-**>(t. t) shown in Fig. 2 is about 10^s, and almost 30 
distinct paths occur. It is thus evident that the occurrence of 
ISI is inevitable at the detector in the case of BU. Furthermore, 
the strengths of the different and independently fading paths 
are in accordance with the scattering function defined in [27] 
for BU radio channels. 

At the base station, the estimates 

d<*> = (d[ k \ 4 k) . - . «?*>) r . e v. v c c, 

k = n = l.-.-.JV. K.N EN (9) 

of the data symbol sequences d (A) defined by (I) must be 
determined. 

Setting out from Fig. 1 , a discrete-time model of the uplink 
in a synchronous JD-CDMA mobile radio system with CRAD 
can be derived in analogy to the single antenna receiver case 
[3], [4], [II], [12]. Each of the A" * K n resulting discrete-time 
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Fig. 2, Absolute value \H a ^° '( r./)| of the simulated complex impulse 
response versus r and t for arbitrarily chosen k and k, t in the case of BU 
and a vehicle speed of 50 km/h. 



radio channels is represented by its channel impulse response 

€t fc = L---.ff.fca = L"-.A" a . 
w = L-* -.IF. K.K a AV eN (10) 

consisting of W* complex samples f$' k " 1 taken at the chip rate 
1/T t .. In the discrete-time model, the data symbol sequences 
d {k) given in (1) are transmitted over A' - K a discrete-time 
channels with the combined channeLimpulse responses 

f A ' a) eC. fc = L---.A\ fc„ = L-».AV 
I = 1..-..Q + W- 1. if. A' a .Q-W €N (ID 

consisting of the discrete-time convolution of h iK ' ka) intro- 
duced in (10) with the signature sequences c w defined by (3). 
ISI arise for W greater than one and MAI occur for IV greater 
than one and/or for nonorthogonai signature sequences c {k) . 
The discrete-time channels with the combined channel impulse 
responses b<*-*« >, Jk = L - - . A\ k a - 1. * - - . K a . according 
to ( 1 1 ) are called (Q + W-l )-path channels. In what follows, 
it is assumed that the combined channel impulse responses 
fc(*.*a) introduced in (1 1) are known at the receiver which can 
be guaranteed by using perfect channel estimation. Channel 
estimation shall not be considered in what follows because 
the effect of channel estimation errors on the data detection is 
not of interest when investigating the potential of different data 
detectors. The reader is referred to [20], instead, where channel 
estimation algorithms applicable to the considered JD-CDMA 
mobile radio system with CRAD as well as their effect on the 
data detection are studied in detail. 

Fig. 3 shows a schematic of the uplink in a JD-CDMA 
mobile radio system with CRAD. The (Q + W - l)-path 
channels are modeled as finite impulse response (FIR) filters. 
Fa^h FIR filter consists of a shift register with (Q + W - 2) 
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Ftg. 3. Discrete-time model of the uplink of a JD-CDMA mobile radio ^ 
system with CRAD. 

Firstly, the well-known case of a single receiver antenna 
with label k tt is considered [4]. According to Fig. 3. the 
received sequence e ( *«> of length {N - Q + W - 1) prevails 
at antenna fc«. Equivalent^ to the single antenna receiver 
discussed in [4], each received sequence e ( * o) consists of a 
sum of A' sequences, each of length (N • Q + W - 1) and 
containing the data symbol sequences d (fc) introduced in (1), 
which are perturbed by an additive stationary noise sequence 

ng^eC. k a = l.---.A a . 
Vi = L---.JV-Q + W-l. ka,N.Q.WeN (12) 

with zero mean and covariance matrix [24], [30] 

R <M(*«) = £{ n ^^) ,r }. fc a = 1. • * * . K tt . K a e N. 

(13) 

Introducing the data vector ^ 
d = ( d d)r -d (2)r . . . d (K) T) r = {dlmd2 ... d K .xf- 

K.N eN H4) 
where the components of d are given by 

K.NtN (15) 

and defining the (N ■ Q + W - 1) x K - N-matrix 

i=l.--.NQ + W-l. 
A (*.) = {A] k ; ) ).j = 1- ■■K -N, (16a) 

k„ = I.- - ■ .K„. 



(»-U+f..V-(A— l)+n 



r &«*■-*■- » for fc= 1. 

fc„ = 1. - • - . K a . 

1 = l.--.Q + W -I. 
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Jfie received sequence can be represented by 



9 



e .v Q+n 



= A^ ) d + n l 



fca = 1." 



.K a . K a .N.Q.W<=\ 



(17) 



Now, the unified mathematical representation for K a re- 
ceiver antennas is presented. With the K a - (N - Q + M — 
1) x K ■ iV-matrix 

A = (A Lj ) = (A^ 1 ". A< 2)r --- A< A '« )T ) r K a eN (18) 
and with the combined noise vector 

n = (n( 1 > r .n( 2 > T ---n^-> I ") r 

= (ni.n 2 "-rtjfc«-(A--Q+ir-i)) • 

K a .N.Q.W€N (19) 



where 

S = (S,J. i.j = l. (25) 
is a square matrix with K • JV rows and A" - N columns, 
M={M tmJ ). / = !-■ 



where 

4l f „(**> 

ri {N-Q-\-W~\Hk 0 -l)+n— Tl n ■ 

n = l.---.iV-Q + VP - 1. 



Ar a — 1 . * • • . k a - 

k a .N.Q.\ve? 



(20) 



bolds, having the covariance matrix 

m R n = E{nn* r } 
£1 / R{, 1)(1) 



x.(2)(l) 
■tin 



VR (K.)(1) 



-tin 



(2)(2) 



R 



t,(A'«)(2) 



= J5?{n< 4 >n<'>* T }, 
i.j = l.^-.A^, AT a eN 



(2 la) 
(21b) 



combined received vector is given by 



jiAd + n. AT^iV.Q.WeN (22) 



with 



d *. { Jk a ) i. _ i ... 1. 



(23) 



./ = L-.-.A' a .(.v-Q + ir - 1) 
is a A' - A r x A" H - (N - Q + W - l}-matrix and 
d = (d\. (I2 •• • &k x) T 



(26) 



(27) 



is the estimate of d defined by ( 14). The choice of the matrices 
M and S determines the equalizer type. The estimated data 
symbols d n contained in d of (27) are either continuous 
valued, referred to as d c , n contained in d c , or discrete-valued, 
referred to as d q , n contained in d q . The continuous valued 
estimates d c .„ must be quantized [4] in order to yield the 
desired discrete-valued estimates </« fc) introduced in Section 
II. Furthermore, it is assumed that (d^> r .d^ ,r ---d (A ' )r ) r 
introduced in Section II is always identical with the discrete- 
valued estimate d q . 

In the following, the representations of M and S in the 
case of CRAD are introduced for the DMF, the DMF-BDFE. 
the ZF-BLE, the ZF-BDFE, the MMSE-BLE. and the 
MMSE-BDFE. It shall be shown for the latter five JD 
techniques that M always contain a DMF. The transmission 
of data symbols d\t ] with E{d { n } } equal to zero is assumed 
in what follows. In this case the SNR 



7(fc.n) 



variance of the sig nal component 

variance of the noise and interference components 



noiseless. 1 isolated 



E{ ^'-(ifc-lH-n 



(28) 



The combined received vector e of (22) is processed in a data 
detector in order to determine the estimates d {k) defined by 
(9). 

Hi. Joint Detection Techniques for 
Coherent Receiver Antenna Diversity 



associated with the estimate d. Y . (fc _ 1)+n of the data symbol 
d { n ] transmitted by user k at time nT s at the output of the 
considered equalizers is given. 

B. Decoirelating Matched Filter 

Conventionally, the data detection in a JD-CDMA mobile 
radio system for a single antenna receiver is accomplished by 
using a DMF followed by symbol-rate samplers, cf. Section 
I. The extension of the DMF to CRAD is straightforward by 
using the definitions for A, R„ and e introduced in Section II. 
With the notation Diag(X ; . £ ) [32] denoting a diagonal matrix 
containing only the diagonal elements of the matrix X and 
with the Cholesky decomposition [31] 

- 1 - L* r L 



r; 



(29) 



A. Bask Concept 

In what follows, it is assumed that the received vector e 
defined by (22) is completely known at the receiver before 
the data detection is carried out. The basic concept of the 
JD techniques considered in this paper is given by the set of where L is an upper triangular matrix 

equations L = (L Lj ). L;j = 0 V / > j. 

Sd = Me (24) /. j = 1. - • - K tt • (.V • Q + W - 1) (30) 
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the continuous-valued output of the DMF for CRAD is given 
by 

d< = A^R^e (31a) 
= {LA)* r Le. (31b) 
= Diuf>([(LA)* r LA],,)d + 



u^efui part 



j(LA)* r LA - Dia g([(LA) ,r LAl,,,)]d+ (31c) 

ISI and MAI 

{LA) +r Ln 

<■ v ' 

perturbation 

where [X]; ^ designates the element in the /th row and the 
jth column of the matrix X, According to (3 1 ), the estimate 
d r contains ISI and MAI as well as a . perturbation term 
associated with filtered noise. It follows from (31b) that S 
is the K - N x A" * iV-identity matrix and M is represented by 



M= (LA)* r (L) 

matched filter decorrelating filter 



(32) 



The operation Le decorrelates, i.e. prewhitens, the noise. 
Therefore, L is a prewhitening or decorrelating filter (DF). 
The signal Le is fed into the filter (LA)* T which is matched 
to the concatenation of the A' - K a discrete- time channels 
with the impulse responses fo^-M of (1 1) and the DF L. The 
concatenation of L with (LA)* r is a DMF. 

For convenience, the Hermitian matrix 



E = A* T R~ 1 A = (LA)* r LA 



(33) 



is introduced. Under the assumption that the data symbols d n 
are samples of a stationary process with covariance matrix [30] 



R (/ = E{dd*< } 



(34) 



the SNR ^{k. n) at the output of the DMF introduced in (28) 
is given by the equation shown at the bottom of the page (35). 

The SNR *)(k.v) introduced in (35) is maximal when 
neither ISI nor MAI are present. In this case, (35) reduces to 

>>[k.n) = E{|ef} * [E] x ,,_ 1K/) . x ,,. 1)+/! . 

fr=l.---.A'.K = l.-.-.AT. (36) 

The conventional JD technique represented by the DMF 
is only suitable in the case where both ISI and MAI are 
negligible, as it is the case in spread spectrum multiple 
access (SSMA) radio systems [1] with low spectral efficiency. 
However, both ISI and MAI are usually considerable in mobile 
radio systems thus leading to impairments. In the follow- 
ing sections, JD techniques with considerable performance 
improvements over the DMF are presented. 



C. Decorrelating Marched Filter Block 
Decision Feedback Equalizer 

In this section, a decision feedback version of the DMF for^ 
CRAD termed DMF-BDFE for CRAD shall be derived setting \ 
out from the DMF for CRAD discussed in the preceding 
section. The DMF-BDFE for CRAD is based on interference 
cancellation, cf. e.g. [10], [22J. Its basic principle is presented 
in the following. For convenience, (31a) is displayed using the 
elements of the vectors d c and e as well as of the matrix M 
of (32) in the following way: 

d CM = M »-J e J' w = l----.*'.-V. (37) 

According to (37), the estimation of d c can be carried out 
recursively with n being the label of the recursion steps. In the 
following, it is assumed that the recursion starts with n equal 
to K • N and ends with n equal to one. Therefore, rfr-.Avv is 
determined during the first recursion step with label n equal 
to A'- N. Based on d c .K \'~ the discrete-valued version d q .K-x 
is generated by quantizing <J c .a'-.v- The effect of those parts 
contained in the received vector e which are dependent on the ^ 
data symbol dKX can now be reduced because the estimate 
dq.K-x of dK-x is known. Hence, d q _K-x is modulated at the 
receiver, i.e. 

is generated, and then (38) is subtracted from the jth compo- 
nent e J9 j = 1. - * • . K a * (N • Q + W - 1), of the received 
vector e at the DMF-BDFE input. In the next recursion step 
with label n equal to (K - N - l)^d q . K x-i is determined. 

Now, the effect of dx-x and (Ik x-i on the received vector 
e is reduced by subtracting the modulated versions 
k x 

^ Ajjd tI , t — Aj.K-X-ldq.K-X-l + Aj.k xdq.K-x- 



i~KX~l 



j = l.---.Jf a -(JV-Q + W-l) (39) 



of dq.K x~\ and d q .K X-i from the jth component e,, j = 
1. • ■ ■ . K a • (N * Q + W - 1), of e at the DMF-BDFE input. 
The described procedure is continued until d q is entirely ^ 
determined. Assuming that a is decremented from A' • .V down 
to one. the mathematical representation of this procedure is 
given b> 

AV(.Y<?+IV-1) / kx > 

dc,n = ]T M„.j \<ij~ A j-> (I <i-> 



(40) 



Introducing the notation Tri(X) for the strict upper trian- 
gular matrix with vanishing diagonal elements consisting of 



E{\<e } \-}-([E\^) 2 



[EIWEj„ + (1 - 2Ro{[ER,/]„. v }) • [E]„.„ + £{!<tf'| 2 } ■ ((£],,„)- 



v = N-{k-l) + n. Jfc = !.•••. A. w = l.-- 



(35) 
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■Jt upper triangular part of the matrix X and introducing the 

//definitions 



ff 



\A\ = (-r.-J-i+i ■■-■r J ) T i<j- 



(41a) 



TO 



+ 1 



(40) can be displayed as 

[Tri{E)]i:"+ 1 A , A .)(([d t .], 1 ,) r . ([d 9 ]^-) r ) r 



= A R 



Tn-l. 



(42) 



with I denoting the K ■ N x K • N-identity-matrix. The matrix 
[Tri(E)]^-"y ^-. v represents the feedback operator. With 

M = A* r R- 1 . (43a) 
S = (M^a-,,- [Tri(E}]^,. v ). (43b) 
•-- -»+^^r (43c) 



a = (([d c ji) r .([d,]E«) r ) a 



NMe structure of (42) is recognized. Furthermore, the 
J0MF-BDFE obviously contains a DMF, cf. (43a). However, 
Itfeth S and d given by (43b) and (43c), respectively, must be 
updated every time instant nT s with n decreasing from K * N 
down to one. The signal processing expense of the presented 
03MF-BDFE is only moderately increased compared with the 
43MF presented in the preceding section. 
~ With the definition 



F = f E-Tri(E). 



(44) 



\ ffie SNR -y{k. n) at the output of the DMF-BDFE is given by 
Mthe equation at the bottom of the page assuming that all the 
^past decisions which are fed back are correct. The SNR ^{k.n) 
Hsf (45) is generally lower than *)(k.n) at the output of the 
DMF given by (36). Error propagation may occur in the case 
of incorrect past decisions. Nevertheless, the performance of 
^ the DMF-BDFE can be improved by applying channel sorting 
similar to [3], (53). 

D. Zero Forcing Block Linear Equalizer 

As mentioned in Section I, a straightforward JD technique is 
represented by the 2F-BLE [3], [4], [24] which is based on the 
Gauss-Markov estimation [23]. The ZF-BLE was discussed 
in [3], [4], 1 24] for the case of a single antenna receiver. 
In this section, the extension of the ZF-BLE to CRAD is 
demonstrated by considering A of (18). n according to (19) 
and e defined by (22). 



In analogy to the single antenna receiver case [3j. [4], [24], 
the ZF-BLE for CRAD minimizes the quadratic form 



Q(d,) = (e - Ad^R-^e - Ad r ) 



(46) 



where d r is a data vector with continuous valued elements 
d c n , n— 1. - - - . A" - A r . The minimum of t^(d r ) is associated 
with the continuous valued and unbiased estimate 



d f =(A' T R; 1 A)- 1 A* r R- 1 e 

= d + (A^R-'Aj^A^R- 1 !! 



(47) 



perturbation 



of d given by (14). According to (47), the estimate d r is 
free of ISI and MAI but still contains a perturbation term 
representing filtered noise. 

It follows from (47) that S is the A' * A r x A" • iV-identity 
matrix and M can be displayed as 



M= (A^Rj^A) 1 A* r R~ 1 



(48) 



In analogy to [3J, M can be further developed. With the 
Cholesky decomposition [31 J 



A* r R~ 1 A = H* r E 2 H 



(49) 



where H is a unit upper triangular matrix and E is a diagonal 
matrix 

H = (H KJ ). H tJ = 0 V / > j. H ut -IV/. 
'/\j = 1. A r . (50a) 

£ = Diag<tr^) (r ( ,6i i = 1- - • - . K ■ A r . (50b) 

(48) can be displayed as 



M = 



(EH)" 1 



(H^Sr 1 (LA) 



ISI and MAI canceller whitening filter matched filtei 

x 00 ■ (51) 

decor relating filter 

As mentioned in Section II1-A, the ZF-BLE for CRAD 
contains a DMF which is applied to the received vector e 
of (22). The output of the DMF is fed into the whiten- 
ing filter (WF) (H* T £)~\ The combination of the filters 
L, (LA)* r , and (ff^E)" 1 shall be termed decollating 
whitened matched filter (DWMF). At the output of the DWMF. 
a maximum likelihood sequence estimator (MLSE) can be 
applied. However, the expense of such a MLSE is prohibitively 
high. Therefore, using a linear ISI and MAI canceller such as 
(EH)"" 1 is more favorable. 



l( Ln)= E{\&\*Um^ 



[FRrfFT]„ .„ - (1 - 2Re{[FR*U}) - [EU + £{|<4, fr) |-} • «EU) 2 ' 
N-{k-l) + n. k=l.--.K. n=l.--.N (45) 
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With E defined by (33), the SNR 7(6. w) at the output of 
the ZF-BLE can be displayed as 



The matrix [H - I]^" v k x represents the feedback opera, 
With 



7(i-») = 



[E 1 ].v.(fr_i) +n ..v.(fc-i) +H 
fc=l.---.A'. n=L*--.iV. K.NeN (52) 

which is in general lower than y(Ln) of (36) at the output 
of the DMR However, the ZF-BLE performs better than the 
DMF when ISI and MAI are considerable. 

£. Zero Forcing Block Decision Feedback Equalizer 

In this section, the extension of the ZF-BDFE [3] to CRAD 
is demonstrated. The ZF-BDFE for CRAD shall be derived 
setting out from the ZF-BLE for CRAD discussed in the 
preceding section. The structure of the unit upper triangular 
matrix H can be exploited for establishing the ZF-BDFE. 
With (47) and (51), the equality 

Hd c = d c + (H - I)d c = IT 1 (H* r S)~ l (LA)* r Le (53a) 
= Me (53b) 

holds, where the matrix 

M = S- 2 (H %r S)~ i (LA)* r L (54) 

is used and I denotes the K N x K - AT-identity-matrix. 
According to (54), the ZF-BDFE contains a DWMR From 
(53b), the identities 

AV(A'-Q+ir-l) 

dc* -at = Mk-x.jZj. (55a) 

K-S (-Y-Q+H--1) 
j=n+l j=l 

n= h--.K-N- 1 (55b) 

follow. According to (55a), the estimate d c h x is deter- 
mined by the linear superposition of the K a - (N • Q + 
W - 1) weighted elements Mk-x.j^j- All the other esti- 
mates rf c .„, are moreover influenced by the weighted es- 
timates H n .„+ 1 d CM +i, H, imn +2tlcM+2"~ Hn*K'Xd c ,K-X- c f- 

(55b). Under the assumption that the set of equations given by 
(55a) and (55b) is solved recursively with n descending from 
K • N down to one, a ZF-BDFE can be realized. Replacing 
d Cm j* j = (n+ 1) • - (K-N). in (55b) by the quantized versions 
d q .j introduced in Section III- A yields the ZF-BDFE with 

d c> K -X — Mk ' N - J ej ' 

(56a) 

KX K a {X-Q+\Y-l) 

ft = 1. - ■ - . A" • JV — 1. (56b) 
With (4la) and (41b), (53b) results in 

(fflk-lv.- (H - ll^.A-)((M) r - iW^ff = Me. 

(57) 



S=([I]^.[H-I]^ X ). (58a) x 
d-(([d,]j/.([d ( K) r ) r (58b) ^ 

the structure of (24) is recognized. Both S and d given by 
(58a) and (58b), respectively, must be updated every time 
instant nT^ cf. Section III-C. 
The SNR y(L n) at the output of the ZF-BDFE is given by 

7(*-«) = E{|f4 fc> | 2 } * (^V-(fr-I)+n„V.(A-l)+«)" 

k = I. - .K. n = l.---.iv T . R\NeN (59) 

assuming that all the past decisions which are fed back are 
correct. The SNR 7(fc.n) of (59) is generally greater than 
7(fc. 71) at the output of the ZF-BLE given by (52). However, 
in the case of incorrect past decisions, the performance of 
the ZF-BDFE suffers from error propagation. As already 
mentioned in Section III-C, the performance of the ZF-BDFE 
can also be improved by applying channel sorting similar to 
[3], (53). * 

F. Minimum Mean Square Error Block Linear Equalizer 

In [3], the MMSE-BLE was investigated for the single 
antenna receiver case. In this section, the extension of the 
MMSE-BLE to CRAD is presented. Again, the definitions for 
A given in (18), n according to (19) and e introduced in (22) 
are used. 

In the case of the MMSE-BLE for CRAD, the quadratic 
form 

Q(d c ) - E{(d c . - d)* r (de - d) } (60a) 
-£{||d c -d|| 2 } (60b) 

must be minimized, cf. [3], [23], [24] for the single antenna 
receiver. The data vector d c , consists of the continuous- valued 
elements d c . n < n = L--.K • N. The quadratic form Q{d c ) 
assumes its minimum when d 0 is equal to the continuous- 
valued and unbiased estimate \ 

d c = (A^R^A + Rj^^A^R^e (61a) 

^(i+fR^A^R-^)" 1 )^ 1 
> v ' 

Wiener filter W 0 

x (A^R^A^A* 7 !*.-^ (61b) 
> ^- ' 

ZF-BLE 

= Diag({W 0 } M )d + [Wo - Diag({Wo}u)]d 



useful part 



ISI and MAI 



+ W 0 (A* T R- 1 A) V r lCn. 



(61c) 



perturbation 



where I designates the K - N x A' - iV-identity-matrix. 
According to (61), the estimate d c contains a useful term 
an ISI and MAI term as well as a perturbation term, and it u 
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? {KT} 



[W 0 B^„.„ - 2Re{[W 0 R << ] l ,„ • [Wo]:.,} + ^{l^l 2 } • |[W 0 ] fAI ,| 2 
i/ = iV-(fc-l) + n. fc = l.---.AT. 
n= iV. K.N €N 



(64) 



the output of the ZF-BLE followed by the Wiener filter [24, 
p. 373] 

Wo = (i + (RrfA* T R~ 1 A) -1 ) ~\ (62) 

Due to (61), S is identical with the K -N x K - JV-identity 
matrix and M can be displayed as 

M = Wo(EH)- 1 (H* T S) -1 (LA)* T L. ((63)} 

According to (63), the MMSE-BLE for^CRAD contmns^ a 
DWMF whicBTTs^applied to the received vector e of (22). 
The output of the DWMEJkXed^ 

cMceller X?H) ~ 1 ~ which, is followed by .the ^ Wiener filter t 
Wo, Since the Wiener filter minimizes the expectation of the 
squared norm of the estimation error vector (d c — d), the 
MMSE-BLE leads to a better performance than the ZF-BLE, 
cf. [3] for the single antenna receiver. Furthermore, it can 
easily be shown that the estimation errors (d Ctft — d n ) and the 
estimated data symbols d c ri are uncorrected at the output of 
the MMSE-BLE for CRAD. 

The SNR -y{L n) at the output of the MMSE-BLE is given 
by the equation at the top of the page (64). which is in general 
greater than j(k. n) at the output of the ZF-BLE introduced 
in (52). 

G. Minimum Mean Square Error Block 
Decision Feedback Equalizer 

In this section, the extension of the MMSE-BDFE [3] to 
CRAD is demonstrated. The MMSE-BDFE for CRAD shall 
be derived in a similar way to the ZF-BDFE by setting out 
from the MMSE-BLE for CRAD discussed in the preceding 
section. With the Cholesky decomposition [31] 

A* r R~ 1 A + Rd 1 = (£'HT r 2'H' (65) 

where the matrices 



Lj = 1.---. K -N. (66a) 
£' = Diag(<7^).<^ € H- * = L • • • . K- N (66b) 

are used, and with the matrix 

M = E'- 1 (H.'* T l?y l (LA)* T L (67) 

the MMSE-BDFE can be represented as follows: 

(ffl^v.„-[H' -i]iS l A -^)((M) r - (K]^-) r ) r 

- Me. (68) 

The matrix [H -IJ^-"vk v represents the feedback operator. 
With 

S=([Ilk- 1 J v.„-[H'-I]^i-..v). (69a) 

d=(([dcli) r ([d,]'A^) r ) (6%) 

the structure of (24) is recognized. 

The SNR 7(Ar. n) at the output of the MMSE-BDFE is given 
by the equation at the bottom of the page (70). assuming that 
all the past decisions which are fed back are correct. The SNR 
7(fe.n) of (70) is generally greater than -y(k.n) at the output 
of the ZF-BDFE given by (59). 

IV. Simulation Results 

In this section, the average uncoded bit error probabilities 
F e obtained with the six JD techniques introduced in Section 
III are presented for K a equal to one and K a equal to two 
receiver antennas in the case of transmission over rural area, 
typical urban and bad urban mobile radio channels. In the 
simulations, a user bandwidth of 2 MHz is chosen and the 
following parameters are used: 

K* = 8, (71a) 
# = 20, (71b) 



7(fc, n) = 



^{Kf}-!^) 2 -!^)" 1 ]-! 



v = N-{k-l) + n. * = !.•■-. A". n = l.---.N. 



K.N 6 N (70) 
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Q = 14. (71c) 

W = 2 (rural area). (7 Id) 

W - 11 (typical urban), (7le) 

W = 21 (bad urban), (7 if) 

r, = ().5mus. (7ig) 



The modulation scheme is 4PSK. The maximal data rate l/% 
is approximately 143 kbit/s. The randomly generated signature 
sequences c^ k \ k — 1. • ■ - , 8, are binary. 

Variations of the received power due to path losses and 
iognormal fading are assumed to be eliminated by applying a 
perfect power control. However, the variations of the received 
power due to Rayleigh fading are present in the received 
signals. In the simulations, 

R<* — I (72) 

corresponding to 

e||4 A) | 2 |-1 (73) 

is chosen. The additive noise is assumed to be white and 
normally distributed with covariance matrix 

R„ = a 2 -I. (74) 
Therefore, the SNR per bit 

£k |lb<"°>|| 2 

Tfc(fc.n) = Tfc(fe)= J] a 2 - 
k a =i 

k = l.. -.h\ n = l. K.K a ,N € N 

(75) 

at the input of the data detector is independent of n. With 
the SNR per bit -y{k,n) at the output of the data detector 
operating with a particular JD technique presented in Section 
III, with the well-known bit error probability Pt{l(k.n)) for 
4PSK for signaling over an additive Gaussian noise channel as 
a function of the SNR per bit 7(fc. n) [2] and with the density 
function p{i{k,n)) of the SNR per bit j{k,n) [2] resulting 
from the varying combined channel impulse responses 
the average uncoded bit error probability is given by 

-rt7(k«))d(7(*.»)) ( 76 > 

[4]. In the case of a fixed choice of b^ ka) and cr 2 , P b (i\k, n)) 
can be determined. However, the variations of b^ ka) must 
be simulated, thus allowing the calculation of P e according 
to (76) [4]. 

The authors intend to determine lower bounds on the error 
rate performance and thus assume that all the K K a combined 
channel impulse responses b (fc,fcft) are linearly independent 
Cuncorrelated"). In practice the K-K a radio channels can, for 
instance, be made uncorrected by choosing sufficiently large 
spacings of the receiver antennas, cf. [33, ch. 9] for a detailed 
discussion. However, even close proximity of the receiver 
antennas leads to uncorrelated radio channels [17] resulting 



from antenna coupling effects. Nevertheless, the presenV 
concept of JD with CRAD is not restricted to the case of tvi\ 
receiver antennas, which are widely spread, but can also be x 
applied to antenna configurations consisting of many closely \ 
spaced antenna elements for which case correlations between 
the received signals occur. Detailed investigations concerning 
the effects of antenna spacing and antenna configuration on 
the correlation properties of the received signals, and thus on 
the data detection by using JD techniques, shall be carried out 
in the future. 

In Fig. 4-6, P e achievable with the six JD techniques of 
Sections III-B-F are shown versus the average SNR per bit 
and per antenna i b /K a obtained by averaging over the SNR's 
76(fc) of (75) under the assumptions that 

• a 2 and b<*- A -\ k = I,- -, K y k a = 1, are 
perfectly known at the receiver and 

• no error propagation occurs in the case of the 
DMF-BDFE, the ZF-BDFE and the MMSE-BDFE. 

In the case of error propagation, which shall be studied in the 
near future, a moderate SNR degradation is expected. 

It follows from Fig. 4-6 that for the considered channel 
models and for one, as well as two receiver antennas, the 
MMSE-BDFE performs best followed by the ZF-BDFE, the 
MMSE-BLE and the ZF-BLE. The average uncoded bit error 
probabilities P e associated with the DMF for -y{k,n) of (36) 
cannot be reached by the MMSE-BDFE, the ZF-BDFE, 
the MMSE-BLE, and the ZF-BLE. Therefore, the DMF 
for 7(fc,7i) of (36) represents a lower bound to P e of the 
considered suboptimum JD techniques. The average uncoded 
bit error probability P e of the MMSE-BLE always merges into 
P e of the DMF for % /K a -* 0 and into P e of the ZF-BLE for 
j b /K a -» oc, cf. (31a), (47), and (61a). The same statement 
is valid for the decision feedback versions of the ZF and the 
MMSE equalizers. 

As expected, P c obtained with the DMF for 7(6, n) of 
(35) and with the DMF-BDFE for 7(fc,n) of (45) merge 
into irreducible error floors because neither ISI nor MAI are 
fully removed. According to Fig. 4-6, the magnitudes of these 
irreducible error floors are dependent on both the considered 
channel model as well as the number K a of receiver an- 
tennas. The effect of the channel model on the magnitudes 
of the irreducible error floors is rather small whereas The 
magnitudes of the irreducible error floors are significantly 
decreased when K a is increased. Since the performance of 
a DMF and a DMF-BDFE is limited to such irreducible error 
floors, neither DMF nor DMF-BDFE are well suited for JD 
in interference limited scenarios. Nevertheless, an acceptable 
error performance is obtainable when powerful forward error 
correction (FEC) coding is applied in combination with the 
DMF for CRAD or the DMF-BDFE for CRAD in the case of 
K a equal to two receiver antennas. 

The application of two receiver antennas instead of a single 
receiver antenna generally leads to a performance improve- 
ment of «3 dB because the unperturbed parts contained in 
e received over the two receiver antennas are superimposed 
coherently whereas the perturbations contained in e are added 
noncoherently at the data detector. The additional performance 
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Fi$. 4. Average uncoded bit error probabilities P. versus the average SNR per bit and per antenna i/,/A'„ for transmission over rural area radio channels 
when DMF. DMF-BDFE. ZF-BLE. MMSE-BLE. ZF-BDFE. and MMSE-BDFE are used for JD in the case of A"„ = 1.2. 



improvements due to the combating of Rayleigh fading can 
be determined by comparing the average uncoded bit error 
probabilities P e of the DMF for 7(4. n) of (36) achievable for 
a single antenna receiver as well as for CRAD. The following 
additional performance improvements due to the reduction of 
the impairing effects because of Rayleigh fading are deduced 
from Fig. 4-6: 

additional performance improvement % 3.5 dB@P e 

= 1(T 2 (rural area), (77a) 
additional performance improvement ^ 1.5 dB@P e 

= 10 ~ 3 (typical urban), (77b) 
additional performance improvement ^ 1.25 dB@P e 

= HT 3 (bad urban). (77c) 
The reason for the different values of these additional 
performance improvements is the varied inherent diversity 
potential of the considered mobile radio channels. For instance, 
the inherent diversity potential of the bad urban mobile radio 
channel is the largest. Therefore, the achievable additional 
performance improvement is the lowest. The following total 
performance improvements obtained with the DMF for 7 (A;, n) 
of (36) result: 

total performance improvement 6.5 dB@P e 

= 1()" 2 (rural area), (78a) 
total performance improvement « 4.5 dB@P e 

= 1(T 3 (typical urban), (78b) 
total performance improvement « 4.25 dB@P f , 

= 10" 3 (bad urban). (78c) 



It follows from Fig. 4-6 that the total performance improve- 
ments achievable with the MMSE-BDFE, the ZF-BDFE, the 
MMSE-BLE, and the ZF-BLE by applying two receiver 
antennas are larger than the total performance improvements 
of (78) because the SNR degradations of the considered 
suboptimum JD techniques apparent in (52), (59), (64), and 
(70) are reduced by applying CRAD. For the same reason, 
the average uncoded bit error probabilities P e achieved with 
the MMSE-BDFE, the ZF-BDFE, the MMSE-BLE, and the 
ZF-BLE when CRAD is considered do not differ as much as 
in the case of a single antenna receiver. 

As an example, the average uncoded bit error probabilities 
P e achievable with the ZF-BLE and the ZF-BDFE in bad 
urban areas are discussed, cf. Fig. 6. In order to obtain P? 
equal to 10~ 3 , %/K a approximately equal to 12.7 dB is 
required in the case of the ZF-BLE for a single antenna 
receiver. However, the required %/K a is only about 6.3 dB 
in the case of the ZF-BLE when two receiver antennas are 
applied thus leading to a total performance improvement of 
about 6.4 dB. Therefore, the total performance improvement 
of the ZF— BLE exceeds the total performance improvement 
of (78c) by about 2.15 dB. When the ZF-BDFE is used P c 
equal to 10" 3 can already be reached for %/K a approximately 
equal to 11 dB in the case of a single antenna receiver 
whereas %/K a of about 5.6 dB is necessary when two receiver 
antennas are applied. The total performance improvement of 
the ZF-BDFE equal to 5.4 dB therefore exceeds the total 
performance improvement of (78c) by about 1.15 dB. 
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V. Conclusion 



In the present paper, six well-known JD techniques, namely 
the DMF, the DMF-BDFE, the ZF-BLE, the ZF-BDFE, 
the MMSE-BLE, and the MMSE-BDFE. applicable to syn- 
chronous JD-CDMA mobile radio systems have been ex- 
tended to CRAD and their viability has been shown by 
simulations. Since the MMSE equalizers minimize the expec- 
tation of the squared norm of the estimation error vector (d e - 
d} ? they perform better than the ZF equalizers. Furthermore, 
the decision feedback versions of the DMF, the ZF, and 
the MMSE equalizers perform better than the corresponding 
linear versions. In addition, the application of CRAD yields 
considerable performance improvements over the single an- 
tenna receivers especially for mobile radio channels with low 
inherent diversity potentials. 
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